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Abstract. 

The theoretical description of trapped weakly-interacting Bose-Einstein 
condensates is characterized by a large number of seemingly very different 
approaches which have been developed over the course of time by researchers with 
very distinct backgrounds. Newcomers to this field, experimentalists and young 
researchers all face a considerable challenge in navigating through the 'maze' 
of abundant theoretical models, and simple correspondences between existing 
approaches arc not always very transparent. This Tutorial provides a generic 
introduction to such theories, in an attempt to single out common features and 
deficiencies of certain 'classes of approaches' identified by their physical content, 
rather than their particular mathematical implementation. 

This Tutorial is structured in a manner accessible to a non-specialist with 
a good working knowledge of quantum mechanics. Although some familiarity 
with concepts of quantum field theory would be an advantage, key notions such 
as the occupation number representation of second quantization are nonetheless 
briefly reviewed. Following a general introduction, the complexity of models is 
gradually built up, starting from the basic zero-temperature formalism of the 
Gross-Pitaevskii equation. This structure enables readers to probe different 
levels of theoretical developments (mean-field, number-conserving and stochastic) 
according to their particular needs. In addition to its 'training element', we hope 
that this Tutorial will prove useful to active researchers in this field, both in terms 
of the correspondences made between different theoretical models, and as a source 
of reference for existing and developing finite-temperature theoretical models. 
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1. Introduction 

1.1. Bose-Einstein Condensation 

One of the most exciting developments in modern physics has been the increasing 
sophistication of experimental techniques to cool, confine and manipulate atoms with 
optical and magnetic fields. In recent years this has culminated in the achievement 
of quantum degeneracy in dilute ultracold gases of bosons [1, 2, 3, 4, 5] and fermions 
[6, 7], allowing the creation of novel many-body systems with unprecedented control, 
tunability. and versatility [8, 9]. 

Indistinguishable particles with integer (bosons) and half-integer (fermions) spin 
differ in how they occupy quantum states. While no more than one fermion can occupy 
each state (known as the Pauli exclusion principle) , the number of bosons in a state 
is unrestricted. The difference becomes most apparent when one cools down a gas 
to low temperatures T, where the de Broglie wavelength X^b oc 1/ VT, becomes the 
same order or larger than the distance between particles. The system then enters a 
quantum degenerate regime. The most spectacular manifestation of this occurs for 
bosons, which below a critical temperature Tc undergo a phase transition, or Bose- 
Einstein Condensation (BEG), where particles tend to macroscopically occupy a single 
quantum state — the condensate. 

For the purposes of the discussion, it is interesting to consider the number of 
particles in the condensate, iVo, as function of temperature T. We consider an ideal 
gas of total atom number N confined in a harmonic potential well of the form Vcxt = 
'm{cu^x^ + fjJyy^ + i^lz^) /2, which closely approximates the traps typically used in the 
experiments. A straightforward calculation [10, 11] gives Nq/N = 1 — (T/Tc)^, where 
Tc ~ O.QAhuJhoN^^^ with Who = (t^s^j^^z)^/^. One sees that at T = Tc, Nq = 0, but as 
T decreases the condensate fraction Nq/N increases, until at T = 0, A'' = A^o and all of 
the atoms are in the condcnsat(\ The inclusion of interactions and the finite size of the 
system changes matters quantitatively, but this simple picture remains qualitatively 
useful [12, 13]. 

Thus, at all temperatures below T^ (apart from the physically unattainable T = 0) 
a Bose-Einstein condensate co-exists with non-condensed particles, which collectively 
make up a thermal cloud. For temperatures very close to zero the thermal cloud can 
be, to some extent, neglected, leading to a relatively simple description in terms of a 
nonlinear Schrodinger equation, also known as the Gross-Pitaevksii Equation [14, 15]. 
Despite its evident simphcity, this equation nonetheless contains much interesting 
physics and provides a good description of many experiments, with much of the 
early theoretical work in Bose-Einstein condensation focusing on solving the Gross- 
Pitaevskii Equation, an effort that continues to the present day [16]. 

Nevertheless, it is important to remember that experiments actually take place 
at finite temperatures. A thermal cloud is always present, and as one increases 
the temperature towards Tc the influence of this on the system behaviour will 
become progressively more important. In some situations the thermal cloud is 
absolutely central, for instance in the problem of condensate growth, or the heating 
of the gas under strong external perturbations. Future applications of Bose-Einstein 
condensation, such as precision measurements based on matter wave interferometry, 
would also benefit from a good understanding of the behaviour of the system at finite 
temperatures. The effects of finite temperature also become particularly important 
for low dimensional systems, where the condensate exhibits fluctuations in its phase. 



CONTENTS 



5 



and the usual picture of Bose-Einstein condensation needs to be revisited. 

Theoretical models of BEC at finite temperature are therefore of tremendous 
value. In order to be useful, such models should ideally accurately represent the 
important physics in a problem, while remaining feasible to solve either analytically, or 
numerically, using current computing resources. The development of dynamical finite 
temperature models that satisfy these conditions has proven to be a great challenge. 
Despite significant progress, research remains active, and much work remains to be 
done. 

A major problem that confronts newcomers to the field of Bose-Einstein 
condensation at finite temperatures, is the large number of seemingly very different 
approaches which have been applied to the problem [17]. This is exacerbated by 
the fact that the major researchers have arrived from different backgrounds, and 
therefore speak "different languages". The aim of this Tutorial is to introduce 
the basic theoretical tools and approximations employed, with the emphasis on the 
physics on which these approximations are based. This is done systematically, starting 
from a brief discussion of the zero temperature theory. In our presentation, we use 
the simplest possible notation which requires minimal prior knowledge, and we also 
give numerous references to equivalent approaches based on mathematically distinct 
formulations. Although space constraints do not allow us to derive all different 
approaches in detail, we nonetheless provide an overview which should allow the reader 
to arrive at a more general understanding of the field as a whole. 

1.2. Basic Formalism 

We begin by reviewing the formalism which forms the basis of our discussion in this 
Tutorial (a more detailed introduction can be found in [10, 11, 18]). Throughout 
this Tutorial, we shall assume that we are dealing with relatively dilute weakly- 
interacting Bose gases, in the sense that the relevant interactions are binary collisions 
between two atoms. Although initial experiments were hmited to this regime, more 
recent experiments with attractive condensates [19, 20], Feshbach resonances [21] 
and molecular BECs [22, 23, 24] require the inclusion of three-body collisions, which 
generally lead to modified scattering and loss of atoms from the system. This Tutorial 
does not include such effects, and for further information the reader is referred to Ref. 
[25] and references therein. 

Within this approximation, the hamiltonian of a closed system of N atoms can 
be written as a sum of two contributions, one (/iq) arising from single-particle effects, 
and the other {V) arising from binary collisions. In the 'coordinate representation', 
this takes the form 

JV N 

H = J2ho{rk) + ^^V{rk,ri), (1) 

fc=l k,l=l 

where the factor of (1/2) ensures that the interactions between every pair of particles 
is only counted once. Actually, in order to formulate a finite temperature theory for 
ultracold gases, it turns out to be much simpler to formulate the problem in terms of 
a different representation, known as the 'occupation number representation' of second 
quantization. A brief review of this, following closely the discussion of Ref. [18], is 
included below for completeness - readers who are familiar with this should proceed 
directly to Sec. 1.2.2. 
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1.2.1. Occupation Number Representation of Second Quantization: A system 
consisting of particles can be fully described by an iV-body wavefunction 
^'(ri • • • rjv, which obeys the well-known Schrodinger equation 



Directly solving the Schrodinger equation in the usual way would be very complicated 
for a system consisting of many particles, as it would require one to keep track of 
particle statistics at each level by using appropriately symmetrized products of single- 
particle wavefunctions, a procedure which is neither practical nor necessary. Instead 
of dealing directly with the many-body wavefunction, one typically reformulates 
the problem in a different representation which greatly simplifies the mathematical 
handling. 

The main idea in this 'new representation' is to exploit the indistinguishability of 
particles in order to keep track only of the number of particles in each energetically 
accessible state of the system. Instead of the usual expansion of the many-body 
wavefunction in terms of a complete orthonormal time-independent basis of single- 
particle wavefunctions, one now defines a new complete orthonormal basis set 
|ni,n2, • • -rioo); where denotes the number of particles in state i (corresponding 
for example to a state with energy £j); while this basis set is infinite, the fact that the 
system has a (fixed) finite number of N atoms implies that (for bosons) there are at 
most N states populated, with all remaining states j unoccupied (i.e. Uj = 0). 

The A''-body wavefunction is thus expanded into this 'occupation number' basis 
set, via the mapping 



where c(ni • • • rioo) denote appropriate expansion coefficients in the new basis. Such 
coefficients must satisfy the following two properties: (i) they must be appropriately 
normalized, such that the probability of finding the system somewhere in configuration 
space is equal to one; (ii) they must incorporate the particle statistics, i.e. be symmetric 
(for bosons) under the interchange of two quantum numbers, to reflect the underlying 
symmetry of the many-body wavefunction. Note that the problem of summing over 
all sets of quantum numbers consistent with a given set of occupation numbers is in 
fact equivalent to the problem of putting A'' objects into boxes, with rii objects in 
Box 1, objects in Box 2, • • • objects in Box i and so on. The main advantage of 
this formulation is that one does not need to explicitly consider symmetrized products 
of single-particle wavefunctions, as particle statistics are automatically incorporated 
within the new formalism. The basis set [rii ■ ■ ■ n^) is time-independent, and thus 
all system dynamics is incorporated in the evolution of the (normalized) expansion 
coefficients c(ni • • • rioo); the latter can be shown to reduce to equations for each set 
of values of the occupation numbers {ni • • • rioo} - i-e. the problem still remains very 
complicated. 

The essence of this approach actually lies in the following realization: The 
equations for the expansion coefiic-ients c(ni • • -n^) generally describe the likelihood 
of particles moving around within the accessible levels, e.g. one particle changing 
its energy by moving from level I to level i. Mathematically, this can be visualized 
as the result of the destruction of a particle in level / and the simultaneous creation 
of a particle in state i; this analogy is very convenient, as it enables us to use the 



ih—^{ri ■■■rN,t) = M(ri • • • r;v, i) • 



(2) 




(3) 
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well-known single-particle annihilation (a) and creation (a^) operators of quantum 
mechanics [26, 27] in the suitably generalized form 

a; |ni • • • rij • • • n; • • • rioo) = \/ni\ni • • • n, • • • (n/ - 1) • • • rioo) 

al\ni ■ ■ ■ Ui ■ ■ ■ ni ■ ■ ■ Uoo) = V«i + • • • (n, -h 1) • • • • • • noo) • (4) 

These operators can be shown to obey the 'standard' bosonic commutation relations 



ai,al 



= 0. (5) 



The number of atoms occupying a particular level i can be obtained within the 
occupation representation formalism via = {Ni), where TV, — (x\ai is the number 
operator, and (• • •) denotes averaging over all single-particle states. This is evident 
from 

7Vi|ni • • • ni_i(ni)nj+i • • • rioo) = a\a^\ni ■ ■ ■ ni_i(nj)nj+i • • • n<x,) 

= aj [^i\ni ■ ■ ■ ni^i{ni - • • • rioo)] 

= a/ {fii - 1) + ly/ni\ni ■ ■ ■ ni-i{ni)ni+i ■ ■ ■ rioo) 

= ni\ni ■ ■ ■ ni-i{ni)ni+i ■ ■ ■ Hoo) , (6) 
leading in the orthonormal occupation number basis to 

(Ni) = (m • • • nj_i(nj)nj+i • • • n^\ni\ni ■ ■ ■ nj_i(nj)ni+i • • • rioo) 

= rii . (7) 

Singlc-particlc effects, such as the action of the potential or kinetic energy 
operator on a many-particle system, can only change the state of individual atoms, 
e.g. shift an atom from one accessible state (say I) to another one of different energy 
(say i). Such a change of state is clearly described by the combined action of one 
annihilation (o;) and one creation (aj) operator, i.e. by the product aja;. Collisions 
included within our theoretical discussion involve only two atoms, with both atoms 
typically emerging from the collision in different states to the ones they were in prior 
to the collision. Such a collision is thus associated with the annihilation of an atom in 
state I (a;) and creation in state i (d|) and the simultaneous annihilation of an atom 
in state k (afc)and creation in state j (a]), i.e. the product of two creation and two 

annihilation operators {alojakai). 

As a result, the original problem concerning the evolution of the many-body 
wavefunction (Eq. (2)), has now been mapped onto the mathematically equivalent 
form 

ih^^\^)=H\^) (8) 
where |^) is defined by Eq. (3), and the hamiltonian is expressed as 

H = Y,{i\ho\l)alai + \ Y,{i3\V\kl)a\a]akdi . (9) 



2 

ijkl 



Here we have introduced the notation 



{i\ho\l) = J dr,p*{r)ho{r)ipi{r) (10) 



with the sjrmmetrized form of the non-local interaction potential taking the form 

1 
2 



{ij\V\kl) = l\{ij\V\kl) + {ij\V\lk) , (11) 
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where 

{ij\V\kl) = j dr j dr'^*{r)if*{v')V{v-v')ipi{v')^k{v) . (12) 

This symmetrization accounts impUcitly for all particle statistics, and demonstrates 
the power of this alternative mathematical technique. Comparing the hamiltonian 

in the occupation number representation, Eq. (9), with the original 'position-space' 
Hamiltonian, Eq. (1), we note the following: Operators appearing in Eq. (1), e.g. 
ho{rk), are now replaced by complex numbers, which are obtained by taking matrix 
elements {i\h()\l) of the original operators over the single-particle eigenstates (pi{r) 
- and similarly for y(rfe,r;). The operator nature of the system hamiltonian must 

be however maintained, and this is achieved by the single-particle operators (aj-^') 
appearing in the definition of Eq. (9), whose role is to remove or add particles in 
particular states of the occupation number representation. 

One may wonder where the averaging (• • •) over the initial single-particle 
eigenstates comes from: We have argued that the original Schrodinger equation (Eq. 
(2)) is mapped onto an equation for the normalizcni and appropriately symmetrized 
expansion coefficients c(ni • • • rioo) in the occupation number basis, with all time- 
dependence included in these coefficients. These expansion coefficients are precisely 
the link between the original 'first quantized' position representation formulation of the 
problem and the equivalent 'second-quantized' representation. In order to obtain an 
equation for the evolution of these coefficients, starting from the original formulation 
of the problem in terms of a complete basis of single-particle wavefunctions {(pi{r)}, we 
must first 'project out' the original basis by averaging over single-particle eigenstates. 
Upon mapping the A^-body wavefunction onto the occupation number representation 
basis, the expansion coefficients thus depend on matrix elements over the initial 
eigenfunctions - more details on this procedure can be found in [18]. 

The hamiltonian of Eq. (9) thus represents, in the occupation number 
representation of second quantization the basic system hamiltonian corresponding to 
the original 'position representation' hamiltonian of Eq. (1). Rather than working 
explicitly with the individual levels in the occupation number representation, it is 
often more convenient for brevity to construct linear combinations of such operators 

i>{r,t) = ^ai{t)^i{r,t), ¥ {r,t) = ^al{t)^*{r,t), (13) 

i i 

which are summed over the complete set of single-particle quantum numbers. These 

quantities, which are operators in the abstract occupation-number Hilbert space, 
are called 'field operators'. Here ^^{r,t) represents the addition of a particle at 
point r and time t, while ^'(r, t) removes a particle. Hence they are known as 
creation and annihilation operators respectively. For bosons, these operators satisfy 
the commutation relations [*(r, t), *'f(r', t)] = 6{r - r') and [*(r, t), *(r', t)] = 
[^^(r,t),^'t(r',f)] = 0. All properties of the full quantum problem are actually 
contained within these Bose field operators, and the problem reduces to identifying 
suitable techniques for extracting the desired information. 

Using Eq. (13) we can thus show that the transition from the coordinate 
representation to the occupation number representation of the system hamiltonian 
is achieved via the transformations 

N 

^Ao(rfe)^ dA\r,t)ho{r)^{r,t) (14) 

fc=i 
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N 

E 



V{Tk,ri)^ j j drdr'¥{T,t)¥{r',t)V{r-T')i!{T',t)i!{r,t) . 



Our subsequent discussion will be mainly given in terms of these Bose field 

operators, as this allows for more compact expressions. However, in order to deal 
with some of the more subtle issues, we will occasionally expand the field operator in 
a suitable basis via Eq. (13). 



1.2.2. Basic System Hamiltonian In the occupation number representation of second 
quantization, the system Hamiltonian can thus be written in terms of Bose field 
operators as: 

H = yrfr*t(r,t)/io(r)*(r,t) 

+ J drdr'if^{r,t)if1{r',t)V{r-r')if{r',t)if{r,t) , (15) 

where ho = —h^V"^/ (2m) + Vci^i{r. t) is an operator for a single particle in the external 
potential Ve:i^t{^,t) (typically a harmonic or periodic potential), and V'(r — r') is the 
exact two-body interatomic potential. 

This Hamiltonian is the starting point for all theoretical treatments of dilute 
Bose gases, and includes both thermal and quantum fiuctuations. All theories 
appearing in the literature arise from distinct approximations to and mathematical 
approaches for dealing with this Hamiltonian. In this Tutorial, we start from the 
simplest such approach, which corresponds to the lowest order mean field theory, 
effectively representing zero temperature. We then gradually build up the complexity 
of treatments by incorporating additional effects step-by-step. 

For dilute gjises at very low temperature, the usual procedure is to make a contact 
interaction approximation 

V{r-r')=g6{r-r') , (16) 

i.e. to assume that the two atoms undergo perfectly elastic local collisions, like two 
billiard balls. The strength of the interaction is given hy g = 4Trh?a/m, where a is the 
s-wave scattering length for a particular atomic species, which can be determined from 
experiments. This is a somewhat idealized scenario, which can nonetheless be put on 
firm ground by a more careful treatment - the origin and validity of this replacement 
will be further discussed in Sec. 3.2. 

Substitution into the Hamiltonian (15) then gives: 

H= J dri''<{r,t)hoi'{r,t)+^ J dr ^'^■(r, i)«'^(r, t)*(r, t)^(r, t).(17) 

The starting point for the dynamics of ultracold gases is to consider the equation 
governing the evolution of the bosonic field operator ^(r,t), i.e. the so-called 
equation of motion. The second-quantized form of the appropriate equation can be 
analyzed in one of three distinct pictures, known as the 'Schrodinger', 'Interaction' 
and 'Heisenberg' pictures, depending on whether the state-vectors, the operators 
corresponding to system observables, or both are time-dependent [18, 26]. Our 
discussion so far (see Eq. (2)) has been given in the so-called 'Schrodinger' picture, in 
which the state vectors are time-dependent and the operators are time-independent; 
in this picture, the solution to the Schrodinger equation at time t, given the initial 
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solution at to, is given by the unitary transformation |^'s(i)) = e~''^(*~*°^/^|^'s(io)) 
whore H docs not contain any explicit time-dependence, and the subscript 'S' has 
been introduced to denote the Schrodinger picture. In the interaction picture, both 
operators and state vectors depend on time. Finally, in the Heisenberg picture 
the state vectors, which can be constructed from the corresponding Schrodinger 
picture state vectors via |^H(i)) = e''^*/^|^'s(t)), are time-independent and all time- 
dependence is contained in the operators. 

While all three pictures can be used interc;hangeably, most subsequent discussion 
will be given in the Heisenberg picture in which the equation of motion of a general 
operator Oh can be shown to obey 



ih 



dt 



On{t),H 



(18) 



In studying the system dynamics, we will actually be concerned with the equations 
of motion of the Bose field operator ^. For the hamiltonian of Eq. (17), this evolves 
according to the Heisenberg equation of motion (henceforth suppressing the subscript 
'H' for compactness) 



ih 



dt 



*(r,t),Fl =/io*(r,i)+5*^(r,i)*(r,t)*(r,t) . (19) 



This equation contains all the information we can hope to extract from the system. 



1.2.3. Separating Condensate and Non-condensate Contributions: In order to 
extract information from this equation, it is convenient to separate the condensate 
contribution, which corresponds to the macroscopic occupation of a single quantum 
state, from the remaining part of the Bose field operator. 

For example, assuming initially for simplicity that a single-particle basis 
corresponds to a good basis set, one can approximate the expectation value (aja;) = 
^wrii, where is the occupancy of the single particle state. Bose-Einstein condensation 
would then correspond to the macroscopic occupation of one of these states, ipo, such 
that Uq = Nq ^ 1. In general, interactions between atoms induce mixing between 
different single-particle states. As a result, it is often beneficial to describe the system 
in terms of dressed basis states, known as quasiparticle states, which will be discussed 
later. 

Assuming just one (suitably identified) state of the system to be occupied 
macroscopically, it is natural to re-arrange the Bose-field operator into two parts [28] , 

*(r,i) = <^(r,t) + 5(r,t), (20) 

corresponding respectively to a field operator for the condensate, (j>, and one for the 
non-condensed atoms, 5. These could either correspond to thermally-excited atoms, 
quantum-mechanical fiuctuations, or atoms promoted into higher energy states due 
to interactions. Although such a split is mathematically exact, approximations are 
inevitably required in the subsequent analysis and when assigning direct physical 
interpretations to these operators. This split is essentially equivalent to the separation 
of the zero-momentum mode in the usual textbook discussion of Bose-Einstein 
condensation in a homogeneous system. 

The above operators can be re-expressed in a general basis set {foi^i} as 

^{r,t) = ao{t)>fio{^,t), S{r,t) = ^ai(t)^i(r,t). (21) 
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At this point it is worth noting the plethora of different symbols used in the 

literature to denote condensate and non-condensate operators, the most common 
alternatives to those used here being ^ for the condensate, and V^', ■)/', or 5ip for 
the non-condensate. 

1.2.4.. Concept of Symmetry breaking: The condensate operator, ^(r,t), corresponds 
to the annihilation operator for atoms in a state macroscopically occupied by A^o ^ 1 
particles. Since the operator for the number of atoms in the condensate is given by 
Nq = ajao, we find No\No) = No\No), i.e this operator acting upon a state with 
a definite number of No condensate atoms returns an eigenvalue A^o- Now, since 
(SoaJ — aQao)/A'^ol-^o) = (l/-^o)l-^o)7 in the limit of a large number of condensate 
atoms, I/Nq 0, and the condensate single-particle operators can be thought of 
as approximately commuting. One could now make the Bogoliubov replacement 
(often referred to as the Bogoliubov approximation) [29], whereby oq — ^/Nq- The 
operator (j){v,t) appearing in Eq. (20) is thus replaced by a complex number 4'{v^t) = 
^/No(po(r,t), often named the "condensate wavefunction" . In this approximation, the 
field operator is simply decomposed as 

4>ir,t) = <j>ir,t) + d{r,t) , (22) 

i.e. all operator dependence is contained in the fluctuation operator S(r,t). 

This is a somewhat drastic approximation which has the direct consequence 
that the physical state described by such a decomposition does not satisfy the 
same symmetries as the original hamiltonian. In particular, although the system 
hamiltonian is invariant under a gauge transformation in the phase of the Bose field 
operator (as it is expressed in terms of an equal number of annihilation and creation 
operators), the wavefunction (l)(r,t) no longer shares this symmetry; technically one 
says that this symmetry has been 'broken'. Breaking of the U{1) global phase 
symmetry (i.e. fixing the condensate phase) leads to a non-conservation of the total 
number of particles (since these two quantities are canonically conjugate). The 
consequence of violation of particle number conservation is evident, since in this 
approximation one assumes that the addition or removal of a particle in the condensate 
does not affect the state of the system, i.e. Nq ±1 ^ Nq for large atom numbers Nq. 
This approximation is equivalent to the statement that the ensemble average of the 
Bose field operator is well-defined and non-zero, i.e. (^'(r, t)) = (p{r, t) ^ 0. This leads 
directly to {S{r,t)) = 0, a desirable property for a fluctuation operator. Given the 
typical large system size, such a subtle issue can, in first instance, be overlooked, thus 
leading to the simplified mathematical formalism reviewed in Sees. 2-4. This provides 
significant insight into the underlying physics, while still yielding excellent agreement 
with experiments. However, from a fundamental point of view, this is an important 
inconsistency that is revisited in Sees. 5-6. 

We should further remark here that, deflned this way, 4){v,t) is a classical field; 
such an approximation bears close resemblance to the conventional treatment of the 
electric field in the theory of electrodynamics, where the quantum description in terms 
of photons is replac;e(l by a c;lassical field. Although we shall here identify (j){r,t) with 
the condensate, in principle this can also include excitations of the system, as long as 
their occupation 1 and quantum fiuctuations are negligible - we shall return to 
this point in Sec. 6.1. 

Using Eq. (22), we can approximate the atom density n(r, t) = (^^(r, t)^{r, t)) = 
nc(r, t) + h{r, t) into two contributions, namely a condensate density where nc(r, t) = 
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|(^(r,t)p and a non-condensate density n(r,t) = {5'^{v,t)5{v,t)). However, in 
most current experiments, quantum fluctuations are largely overwhelmed by thermal 
fluctuations. For example, even at T = 0, the typical quantum depletion in 3D is given 
by n — no = (8/3\/7r) Vna^ [12], which is typically only a small fraction of the total 
number of atoms, as the 'diluteness parameter' na^ takes typical experimental values 
around na^ ~ 10""^ (where n denotes the density and a the s-wave scattering length 
discussed later). One thus tends to identify (5(r, t) as the operator for the thermal cloud, 
in which case n(r, f) describes the density of the thermal atoms. Note however that 
quantum fluctuations may become important in low-dimensional geometries giving 
rise to novel phenomena that lie beyond the scope of this Tutorial [30] . 



1.2.5. Basic Hamiltonian Contributions: The 'standard' procedure in modelling 
Bose gases theoretically is to separate the Bose field operator into condensate 
and non-condensate parts and then break down the full system hamiltonian into 
various contributions based on the number of condensate and non-condensate factors 
contained in each of them. In particular, substitution of Eq. (22) into the system 
Hamiltonian (17) leads to 

H = Ho + Hi + H2 + H3 + Hi. (23) 

where 

5, 



Hr 



dr 



Hi= dr dHho + 9\cl>\^)cf>+<l>*{ho + 9\^\^)S 



= Jdr [jt (ho + 2g\ct>\'') 6 + ^{{rfSS + ./.^^t^t)' 



H: 



dr 



H.= lldrSmS, 



(24) 
(25) 
(26) 
(27) 
(28) 



where Hq has no 'hat' as there are no operators within it (it is a purely classical 
quantity represented by a complex function). Consideration of these separate 
contributions forms the basis of our subsequent discussion. In particular, we shall 
show how inclusion of different contributions to the Hamiltonian (23-28) can be used 
to derive progressively more sophisticated treatments. 

We shall present various approaches, which can be roughly grouped into three 
different 'classes': (i) perturbative approaches based on symmetry-breaking; (ii) 
perturbative number-conserving approaches; and (iii) unified number-conserving 
approaches which include fluctuations around the mean field. 

Within the c;ontext of mean-field theory treatments, there are basic;ally two 
equivalent ways of formulating the problem: When one is interested in static 
properties, one can either write down the corresponding equation of motion and 
solve it in the time-independent limit, or, equivalently, diagonalize the appropriate 
hamiltonian in the grand canonical ensemble H—fiN, where /x is the chemical potential 
of the system and N the total number operator. Note that as symmetry-breaking 
approaches violate particle number conservation, any subsequent calculations should 
be performed within the grand canonical ensemble [31]. For dynamical properties, one 
either works with the full equations of motion, or after diagonalizing the hamiltonian. 
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one can consider linearized equations arising from the addition of time-dependent 
fluctuations around the cquihbriuni vahics. 

The description of finite temperature Bose gases has historically relied on 
diverse approaches, with some of the key early papers in this field given in Refs. 
[32, 29, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 28] 
(although this list is by no means exhaustive). Existing theoretical models are based on 
similar underlying principles, and the aim of this Tutorial is to focus on the underlying 
physics while developing the theoretical understanding of such systems from first 
principles. While this Tutorial is aimed at the non-expert, active researchers in this 
field may also find it useful, as we have placed particular emphasis on demonstrating 
the relation between the different theoretical formulations currently employed. 

1.3. Overview 

We present here a somewhat detailed overview of this Tutorial: 

Sec. 2 focuses on the simplest possible mean field theory. This is based on 
the approximation H Ki Hq, in which the non-condensate field operator is neglected 
completely. This gives rise to the well-known Gross-Pitaevskii equation, which is valid 
at T = and in the limit of negligible quantum fluctuations. By linearization we then 
derive equations for the collective modes around the ground state. These are shown 
to be identical to the modes of Bogoliubov quasiparticles found by diagonalization of 
the flrst three terms of the Hamiltonian, i.e. H Ri {Hq + Hi + H-z)- 

Consideration of the flnite temperature case necessitates the identiflcation of 
suitable 'generalized' mean flelds presented in Sec. 3. This requires the inclusion of 
contributions from {H^ + H^) to the system Hamiltonian. 'Suitable approximations' 
can be used to reduce these additional hamiltonian contributions to quadratic form, 
thus introducing flnite temperature mean fleld corrections which lead to diverse 
theories in suitably identifled limits. In the so-called Hartree-Fock-Bogoliubov 
approximation, the applications and limitations of which will be discussed and 
assessed, one typically assumes that the Bogoliubov quasiparticle states are thermally 
populated, but their collisions are neglected and the thermal cloud remains in static 
equilibrium. Technical subtleties associated with the inclusion of effective contact 
interactions are also briefly addressed in our formalism. 

Sec. 4 generalizes this treatment to a dynamic thermal cloud. This requires a 
more careful consideration of the contributions (ifs -|- H^), beyond the usual mean 
fleld. Such a description is important, as in most experiments, the condensate and the 
thermal cloud are not in equilibrium, but continuously exchange particles, within the 
constraint of a constant total atom number. This leads to a more realistic description 
in terms of a dissipative Gross-Pitaevskii equation for the condensate coupled to a 
Quantum Boltzmann Equation for the non-condensate; such an approach includes 
collisions between the quasiparticles and thus accurately describes the known damping 
mechanisms. 

Sec. 5 addresses some subtle issues associated with fluctuations in the condensate 
phase and number-conservation. The absence of Bose-Einstein condensation in low- 
dimensional geometries is briefly discussed and an appropriately modifled mean fleld 
theory is highlighted which can describe the equilibrium properties of such systems at 
flnite temperatures. The discussion of Sees. 2-4 is then revisited under the physically- 
relevant constraint of a constant total particle number. 

This Tutorial would not be complete without a brief presentation of alternative 
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more advanced techniques for dealing with both thermal and quantum fluctuations 
given in Sec. 6 (sec also [52]). These include classical field methods (already reviewed 
in [53]) and stochastic methods. Sec. 7 discusses the existing implementation of the 
presented approaches in two key areas of ultracold Bose gases, thus clearly highlighting 
the benefits and limitations of the diverse approaches presented throughout this 
Tutorial. Sec. 8 summarizes and links the different methodologies, with emphasis 
on the current status of the field, and our progress towards a 'complete' description 
of ultracold Bose gases. 

This Tutorial also contains three fully-linked appendices which discuss certain 
lengthy theoretical models, whose physical content and/or simplified limits are dealt 
with in the main body of the Tutorial. For maximum clarity, the end of each 
section features a brief summary of the main results presented in that section, with 
the essential content of the different theoretical models additionally summarized in 
appropriate tables. 

To aid the reader in selecting the most appropriate sections for their study, we 
note that researchers who merely wish to obtain an overview of the underlying physical 
issues without the need for precision and in-depth understanding, should focus their 
attention on Sees. 2-4 which constitute the main part of this Tutorial. 



2. Zero Temperature Mecin Field Theory 

We start by briefly reviewing the zero temperature formalism in order to facilitate 
easier comparison to our subsequent finite temperature discussion. 



2.1. Consideration of Hi^: The Gross-Pitaevskii Equation 



In the T = limit (essentially) all of the particles are in the condensate, so that 
N = Nq and the noncondensate operator can be neglected {S = 6^ = 0); in other 
words, we set ^'(r,t) cf){r,t). Hence, ^'^(r,i) = (f>*{r,t), and the exact Heisenberg 
equation of motion (19) reduces to the so-called Gross-Pitaevskii Equation (GPE) 
[14, 15], 

_d(p{r, t) 



ih- 



dt 



-^V^ + VUr,t)+9mr,t)\' 



<j>{r,t) . 



(29) 



This equation is a nonlinear Schrodinger equation, corresponding to the zero 
temperature hydrodynamic description of Bose gases, first introduced to study 
vortex lines in an imperfect Bose gas [14, 15]. This equation is analogous to 
the equation describing the electric field in Kerr nonlinear media, only in our 
present context the nonlinearity arises from atomic interactions. Moreover, this 
equation is mathematically analogous to a Ginzburg-Landau-type approach [18], valid 
near the critical regime - although the origin of the various contributions and the 
physical interpretation here is actually quite distinct. Remarkably, the GPE, whose 
first numerical implementation to dilute weakly-interacting trapped Bose gases was 
undertaken by Mark Edwards, Keith Burnett and collaborators [54, 55, 56], provides 
a good description of the dynamics of a Bose-Einstein condensate for a large range of 
problems at temperatures as high as T « Tc/2. 
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One can find static solutions by eliminating the 'trivial' time-dependence via the 
replacement: 

<^(r,^)=.^o(r)e-'''*/^ (30) 

where /i is the chemical potential (which is equal to the energy per particle only for 
non-interacting particles). Substituting (30) into (29) yields the time-independent 
GPE: 

-|^V2 + Fext(r)+5l<^o(r)|^ 



M</'o(r) 



(/.o(r). (31) 



In the form given above, the wavefunction is normalized to the total number of 
particles, i.e. / dr|(/)(r, t)p = N . However, various authors work instead with a 
rc-dcfincd wavcfunc;tion, c6(r , t)/^/N, which is normalised to 1, for which the prefactor 
of the nonlinear term in the GPE becomes gN. 

Note that the GPE can also be derived from variational considerations as follows: 
In general, the condensate energy can be written as a function of the condensate 
wavefunction ^(r) , which is termed a 'functional', denoted by E[(j)]. In the limit of a 
large number of atoms, this takes the approximate form 

fi2 _ , 1 



m = Jdv 



2^ |V0(r)|VFe.t(r)|</.(r)|2+-5|0(r) 



(32) 



The desired optimal form for the condensate wavefunction is actually obtained by 
minimizing this energy functional, E[(j)\, with respect to the wavefunction (f>(r); this is 
achieved by the process of functional differentiation, denoted by SE/6(j), whereby one 
examines how the entire functional E[(j}\ changes as a result of small changes in the 
function 0(r). As the wavefunction is a complex quantity, this minimization should 
consider both real and imaginary contributions as independent, i.e. consider changes 
in E[(j)\ arising from independent variations of (/)(r) and (j)*(r) [57]. In performing this 
minimization, one should impose the constraint of a fixed total number of atoms N; 
this can be included by the technique of Lagrange multipliers: in other words, instead 
of minimizing E, one should minimize (or find saddle points of) the quantity (E— fiN), 
where the Lagrange multiplier /x is identified as the chemical potential of the system 
[10]. 

2.1.1. Hydrodynamic Description Since the wavefunction is a complex quantity, 
one can use the Madelung transformation (t>{Y,t) = i/no(r, t)e'^^'''*^, along with the 

identification of v(r,t) = {h/ni)[Vd{r,t)] as the superfiuid velocity to re-express Eq. 
(29) in the form of a 'conservation of mass' (or continuity) equation 

^ + V-(nov) = (33) 
coupled to a 'generalized Euler equation' given by 

I + V . v) V = -V/.0 = -V (- + Fext + .no) . (34) 



In the above equations, the subscript '0' has been used to highhght the fact that 
all the atoms are in the condensate. Eq. (34) depends on the gradient of the T = 
chemical potential, /Uq, and can be re- expressed in terms of a force F = — (l/m)V14xt) 
and the gradient of a 'quantum pressure' contribution Pq; the latter, defined by 
Pq = {l/2)gnQ — (l/4)noV^(lnno), contains two contributions: the first corresponds 
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to the usual pressure term for a gas, while the second term describes the kinetic energy 
contribution arising from spatial variations of the wavefunction occurring within a 
small scale [10]. 



2.2. The Bogoliubov equations 

The GPE, Eq. (29) also admits time-dcpcndcnt solutions, that correspond to collective 
modes or elementary excitations of the system. These can be found using: 

c^ir,t) = e-^'^'/'^[Mr) + mr,t)], (35) 

where 64> represents fluctuations of the condensate wavefunction around the ground 
state (po- If these fluctuations are small, Scj) <^ (()o, then one can neglect powers of Scl) 
greater than one, a procedure known as linearization. Substituting (35) into (29) and 
subtracting (31) then yields: 

ih^^6^iv,t)= [/io + 25|0oI'-m] ^0(r,t) + #^(5</'*(r,t) . (36) 

The linearized equation (36) can then be solved to find collective modes of the 
system. This is usually done by looking for solutions of the form: 

5</.(r,i) = J2 h.(r)e-^"-* + <(r)e-'*] , (37) 

i 

where i labels the different modes, each with frequency LUi. Substituting (37), together 
with its complex conjugate, into (36), and collecting together prefactors of e"**^'* and 
giwit yields the two coupled equations: 

ho + 25|<^o(r)|^ - n] Ui{r) + g[Mr)]''vi{r) = eM^) (38) 

ho + 2g\Mr)\^ - m] v^{r) + g[(j)*{r)]^ u,{r) = - e,v,{r) (39) 

where = hiVi denotes the dressed (quasiparticle) energy of level i; as discussed 
explicitly below, one actually views the excitations above the ground state as dressed 
particles, or quasiparticles, where the dressing arises from interactions via mean field 
coupling. These are known as the (zero-temperature) Bogoliubov equations, and were 
first used by Pitaevskii in 1958 [15] to discuss vortex excitations in a superfluid. 
Similar coupled equations arise in the Bardeen-Cooper-Schrieffer (BCS) theory of 
superconductivity for the description of Fermi quasiparticles [58]; in this context, 
these equations are often referred to as Bogoliubov-de Gennes equations, a terminology 
which is sometimes also used for their bosonic counterparts. 

Eqs. (38)-(39) can also be derived by an alternative, yet equivalent, procedure of 
diagonalizing the quadratic hamiltonian {Hq + Hi + H2), as discussed in Sec. 2.2.3. 

The BogoliTibov equations are often cast in matrix notation as [56] 



(40) 



i(r) M(r) \{u,ir)\^ f u,{v) 
-M*(r) -L*(r) J \ Vi{r) J ' { Vi{r) 

with 

L = ho + 2g\Mr)f-fi (41) 

M(r) = g[Mr)? (42) 

such that L*{r) — L{r). Self-consistent numerical solution of the GPE (Eq. (31)) and 
the Bogoliubov equations (Eq. (40)) fully describes the static properties of the system 
at this level of approximation [56]. 
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2.2.1. Uniform condensates: In order to explain the notion of a bosonic quasiparticle, 
it is instructive to examine tlic solution of Eqs. (31) and (40) for a uniform condensate. 
In this case, Vext = 0, and the density Uco = |</>oP thus becomes independent of 
position, so that Eq. (31) gives /U = grica- Substituting the plane wave solutions 
Ui{r) = Upe'P""/^ and Vi{r) = t;pe'P"'/^ into (40) yields: 



e(p) = hiop 




2m 



2gn 



cO 



(43) 



hio-p introduces 



Here £p = |pp/2m denotes the energy of a free particle and e(p) 
a new dressed energy corresponding to the energy of a quasiparticle excitation. In 
particular, Eq. (43) yields the Bogoliubov dispersion relation [29]. At large momenta 
p > mc (or equivalcntly, e(p) > fi), this excitation energy approaches that of a free 
particle, e(p) = |pp/2m+(7nco = Sp + gn-co- However, in the opposite regime of small 
p, the spectrum approximates a phonon dispersion e(p) — > cp, where c = ^Jgnoijm is 
the speed of sound. 



2.2.2. Non-uniform condensates: Solving Eqs. (31) and (40) numerically for the case 
of trapped condensates [56], leads to a discrete spectrum of excitations, which at low 
energies correspond to bulk oscillations of the condensate. For example, in harmonic 
traps there exists a set of dipole modes, a,t co = tOx = tOy = tOz, which correspond to 
a ccntrc-of-mass oscillation of the condensate within the trap. Other common modes 
involve changes in the condensate widths, either in-phase in the different directions 
("breathing modes") or out-of-phase ("quadrupole modes"). However, for excitations 
whose wavelength is smaller than the length-scale of variations of the condensate, this 
leads to a locally homogeneous condensate with a quasi-continuous spectrum, resulting 
in sound wave propagation, as observed in [59]. A more complete discussion of the 
different collective modes can be found in [12, 60, 61]. As discussed in Sec. 7.1, the 
study of such low-energy modes has provided historically a systematic test for the 
development of a consistent theoretical formalism for finite temperature effects, with 
Eqs. (31) and (40) providing an excellent description at sufficiently low temperatures. 



2.2.3. Diagonalization of {^Hq + Hi + H2j : The Bogohubov equations (40) can also 
be derived from a diagonalization of the first three terms of the Hamiltonian of Eq. 
(23), i.e. upon approximating H w (iJg + Hi + H2). Due to the non-conservation of 
total atom number, this diagonalization should be performed in the grand-canonical 
ensemble, i.e. when the system hamiltonian H is replaced hj K = H — iiN, where 
N = J dr^^'i' is the number operator. We thus repeat the procedure of Eqs. (24)- 
(26), and separate K into its constitutent parts Kq + Ki + K2- Upon noting that the 
wavefunction 4>o(v) is now time independent, i.e. ^'(r,i) = 4>oi'*^) + ^(r,t), we obtain: 

5, 



dr 



(k*o (ho - /u) 



:l<^o| 



K2 



= j dr [^t {ho + 5|(AoP - m) '/'o + (^0 + 5l<^o|^ - m) ^ 



dr 



(44) 
(45) 

(46) 



Minimization of Kq at constant chemical potential is equivalent to the minimization 
of the energy functional performed in Sec. 2.1 under the constraint of a fixed total 
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atom number, and thus leads to the time-independent GPE of Eq. (31). Note that, 
in performing such minimization, wc have already assumed that a unique minimum 
exists in which the condensate chooses a definite (yet random) phase, consistent with 
our assumption of symmetry-breaJcing in Eq. (22). 

How docs the inclusion of further contributions modify this picture? Firstly, 
regarding Ki, substitution of the GPE into Eq. (45), leads to Kx = 0, upon noting that 
fiQ is hermitian. Moreover, the K2 term (46) can be diagonalized with the Bogoliubov 
transformation [18, 69, 29, 28]: 

5{v,t) = J2 [Mr)m + <(r)/3j(i)] , (47) 

where /3 and are annihilation and creation operators for quasiparticles§ which obey 
the Bose commutation relations [/3i,/3j] = Sij and [fii,Pj] = [/3|,/3j] = 0. Expressing 
the non-condcnsatc operators in terms of quasiparticles via Eq. (47) is analogous 
to the expansion around the condensate mean field given by Eq. (37). From the 
commutation relations for the noncondensate field operator [(5(r), (5^(r')] — 6{r — r') 
and [(5(r), (5(r')] = [5^(r), ^^(r')] = 0, one can then derive the orthonormality relations 
(see also Sec. 5.3): 

I dr [<(r)n,(r) - <(r)^;,(r)] = S^j, (48) 

J dr [u^{r)vj{r) - v,{r)uj{r)] = (49) 

Substituting (47) into (46), one can show that if Ui(r) and Vi{r) satisfy the Bogoliubov 
equations then: 



^2 = ^ CiPjpi -Y.'ij dv\vi{v)\\ (50) 



Since the second contribution in Eq. (50) is typically small (assuming negligible 
quantum depletion), the problem essentially reduces to a system of non-interacting 
quasiparticles with an energy spectrum obtained from the Bogoliubov Eq. (40). 



2. 3. Brief Summary 

The T = behaviour of dilute weakly-interacting Bose-Einstein condensates is fully 
described by the hamiltonian Ho, whose purely mean field nature gives rise to a non- 
linear Schrodinger equation known as the Gross-Pitaevskii equoMon. This is given in 
full by Eq. (29) with its corresponding time-independent expression given by Eq. (31). 
Excitations on top of the condensate ground state can either be described by linearizing 
around the condensate wavefunction via Eq. (37), or equivalently by diagonalizing the 
quadratic hamiltonian {Hq + Hi + H2) in the grand canonical ensemble by means of 
the Bogoliubov transformation of Eq. (47). 



Having established our essential notation, we now proceed to discuss how the 
effect of temperature can be introduced into our treatment. A significant part of our 
subsequent discussion will be to examine how the GPE and Bogoliubov equations arc 
modified at finite temperatures. Importantly, while the GPE is modified by additional 
mean field potentials and kinetic contributions, the structure of the Bogoliubov 

§ The quasiparticle operators are often denoted in the Uterature as and d| . 
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equations is essentially maintained, with modifications taking the form of additional 
contributions to the operators L and M. 

3. Finite Temperature Mean Field Theory: Static Case 
3.1. Approximate Inclusion of 

Wc now extend the formalism to finite temperatures, by explicitly retaining the non- 
condensate operator 5{v,t) in Eq. (22). We proceed via the technique of hamiltonian 
diagonalization already discussed in Sec. 2.2.3. In order to go beyond the zero- 
temperature discussion, we must additionally include the remaining contributions 
(ffa + Hi) of Eqs. (27)-(28) to the system hamiltonian. This Section discusses 
their approximate mean field inclusion, which amounts to a static thermal cloud 
approximation. In Sec. 4 we shall see that a more careful perturbative analysis of 
such terms leads to the additional incorporation of the dynamics of the thermal cloud. 

3.1.1. Conventional Mean Field Approximations: In order to reduce the full 
hamiltonian of Eq. (23) to a desired quadratic form, one may consider imposing the 
following mean-field approximation for the non-condensate operators in H^, 

PPUc^A{p5)P5+{PP)U+{U)PP-[2{p5){P6)+{U){PP)].{bl) 

This approximation is motivated from Wick's theorem [18, 69], which states that 
at equilibrium, an average over multiple operators can be approximated by sums of 
averages of pairwise contracted operators, i.e. 

(,5^^t^^) ^ 2(Jt(5) ((5^,5) + {55) {PP) . (52) 

The above approximation maintains correlations of non-condensate operators only 
to quadratic order. One may thus also wish to approximate products of three non- 
condensate operators appearing in by their corresponding 'quadratic forms' 

P55^2{P5)5 + P{55), PP5^2p{P5) + {PP)5. (53) 

Since by construction (^^^^) = 0, the approximation of Eq. (53) implies that {P55) = 0. 

In Sec. 4.4 we shall show that imposing the approximation of Eq. (51) is equivalent 
to ignoring collisions between thermal atoms, while Eq. (53) amounts to ignoring 
particle-exchanging colhsions between condensed and thermal atoms. We have already 
seen that the system may be split into two sub-components, namely the 'condensate' 
contribution, nc(r, i) = |(/)(r,f)p and the 'thermal cloud' n(r,t) = (P {Y,t)5{Y.t)), 
satisfying nTOTAL(r,t) = (*I'^(r, f)^(r, f)) = nc{Y,t) + n{r,t) . Similarly, the 
approximations of Eqs. (51), (53) motivate the definition of an additional mean field 
contribution fh{r,t) = {S{r,t)5{r,t)). This is often referred to as the pair anomalous 
average, and bears its name from the fact that there is an unequal number of creation 
and annihilation operators being averaged over. An analogous correlation plays a 
dominant role in the BCS theory of superconductivity [58] , where fermionic atoms pair 
up to form so-called Cooper pairs. In the case of Bose-Einstein condensation of neutral 
bosonic atoms, where the condensate mean field (j)(r) is the dominant parameter, 
the anomalous average plays a more minor role for effectively repulsive interactions 
between condensate atoms; such a contribution, which does actually become crucial in 
the presence of attractive effective interactions and molecular BECs, should however 



CONTENTS 



20 



not be dismissed a priori even for repulsive interactions, where it will be shown to lead 
to modifications of the atomic scattering. 

With the mean field approximations (51) and (53), the and parts of the 
Hamiltonian can be rewritten in terms linear or quadratic in 5, as 



9 



dr 



SHi 



= g j dr (20(^t^)^t + ^.c.) + g j dr 



5 + h.c.) (54) 

(55) 



5Hq = 5H^' +5Hi°° = 



-g j dr{p5){P8) - | ^ dv{58){PP){m) 



6H2 = SH^^ + SH^""^ 
= 2g f dr{6U)S^^ 



dr 



(57) 



Note that as the following discussion relies a lot on the inclusion of these and other 
subsequently identified related contributions to the quadratic system hamiltonian, 
the definition, origin and interpretation of all such contributions is suininarized in 
Table 1. The inclusion of the above contributions into the system hamiltonian lead 
to modifications to the original hamiltonians Hq, Hi and H2 of Eqs. (24)-(26): 
5Hq merely introduces a shift in the overall system energy, and is therefore often 
neglected, whereas 5H\ and 5H2 introduce crucial modifications to the governing 
system equations. 

In the above expressions h.c. stands for hermitian conjugate, and we have 
separated off the so-called Hartree-Fock (HF) terms involving the collision of one 
condensate and one thermal atom from the Bogoliubov (BOG) terms which correspond 
to the collisional promotion of two condensate atoms to thermal modes (and 
its inverse process). The physical significance of such contributions is displayed 
diagrammatically in Fig. 1 (top). Such a separation enables us to identify various 
distinct approximations. 



3.1.2. The Hartree-Fock (HF) Limit: In the Hartree-Fock limit, one first discards 
from the quadratic hamiltonian (iJg + Hi -\- H2) the latter contribution of Eq. (26) 
which contains two annihilation or two creation operators. One is thus left with the 



reduced zero-temperature hamiltonian [Hq 
the notation 



Hi -\- H2), where we have introduced 



H^ = J drS^ (ho + 2g\<i>{r)f')s. 



(58) 



This hamiltonian is then generalized by the additional inclusion of the contributions 
{6H^^ + 8H^^ + 5H^^) arising from the 'higher-order' hamiltonians H3 and H4. 
Diagonalization of the resulting hamiltonian in the grand canonical ensemble leads 
essentially to 



K 



HF 



drS^ 



ho + 2g {\(t)o\^ + fio) 



(59) 
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SCHEMATIC OF CONDENSATE - THERMAL ATOM COLLISIONS 



HARTREE 

(Direct) 

Propagators: 



FOCK 
(Exchange) 



■■O" CONDENSATE 
_r^^ THERMAL 



BOGOLIUBOV 
(Pair Excitation) 

Exact ^-;;) 
Interaction: " '"' f 



SCHEMATIC OF T-MATRIX 




+ 



Figure 1. (colour online) Top: Schematic of 'Hartree', 'Fock' and 'Bogoliubov' 
coUisional processes involving both condensate and thermal atoms. Bottom; 
Schematic of the definition of the Transition (T) matrix in terms of the exact 
interatomic potential V{r — r'). In the simplest case, the upgrade is to 
the two-body T-matrix, which can then be justifiably approximated by the 
pseudopotential (along with an upper momentum cut-off). However, both the 
intermediate propagators (denoted by filled blue arrows in the latter term) and 
their corresponding energies may be dressed, either due to thermal occupation 
affecting the scattering into these states, or due to the states themselves being 
dressed to quasiparticle ones. 



with the condensate wavefunction satisfying the generahzed time-independent GPE 

ha + g\(f)o\'^ + 2gha (/)o = M0o • (60) 

The absence of any contributions involving two Hke (creation or annihilation) operators 
in the system hamiltonian in this limit thus implies that the system is still described by 
single-particle energies; however, these energies are modified both by the condensate 
mean field \(t>o\'^ and by the thermal atoms tt-q, leading to dressed 'Hartree- Fock' 
energies of the form 

£,(r) =e, + 2g + no(r)] - ^l , (61) 

where is the corresponding single-particle energy in a trap. The equilibrium thermal 
cloud density is given by no(r) — l'<'i('')P(^l^i) where 'Pi(r) are the corresponding 
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Position z (1 ) 



Figure 2. (colour online) Left: Condensate fraction as a function of reduced 
temperature as predicted by Hartree-Fock theory (solid, black) and for the ideal 
trapped Bose gas via A^o /N = 1 — (T/Tc)^ (dashed, red), where Tc is the transition 
temperature for an ideal gas. Right: Total atomic density profile (solid, black) 
consisting of a condensate (dot-dashed, blue - calculated here within the Thomas- 
Fermi approximation [10, 12]) and a thermal (dashed, red) contribution. Shown 
is the column density of a dilute isotropically confined trapped 3D ^"^Na Bose gas 
consisting of 640, 000 atoms, with a condensate fraction of 20%. (Figure created 
by Stuart Cockburn). 



eigenfunctions, {alai) — [e'^^'^'") — 1]""^; and /? — 1/(A:bT) is the inverse temperature 
(fcs is Boltzmann's constant). 

Let us now briefly remark on the origin of the terminology 'Hartree-Fock': The 
symmetrization of the matrix elements performed within the second quantization 
formalism via Eq. (11) implies that in considering collisions involving thermal atoms, 
one is automatically including both 'direct', or Hartree, contributions of the form 
(Oj|y|Oj) and 'exchange', or Fock, terms (Oj|y|jO) (see Fig. 1 (top)). Careful 
consideration (for a detailed discussion see [10]) shows that the appearance of the 
exchange term, which is not present for a pure condensate, implies that for a system 
with a fixed number of atoms, the interaction energy above the critical temperature is 
double that at T = 0. Hence, the relative factor of two in the interaction contributions 
appearing in Eq. (60), whereas such a distinction does not arise in the respective 
equation (Eq. (61)) for the excitation energies above the ground state, where the 
mean fields can be thought of as providing an additional potential for the atoms. 

The Hartree-Fock theory was used in numerous early studies of homogeneous 
condensates [62, 63]. To proceed further with using this theory in trapped gases, 
one often imposes an additional semi-classical approximation [10, 64] which assumes 
that all relevant quantities of the trapped system (such as densities) vary slowly on 
the typical length scale of the confining potential (i.e. the system locally resembles 
a homogeneous gas). The excitation spectrum can then be described in terms of 
momentum p via 

£(r,p)=e(p)-HKxt(r) + 25[|0o(r)|' + no(r)] - fi , (62) 

where £(p) — |p|^/2m is the energy of a free particle of momentum p. This 
semiclassical approximation enables a direct calculation of the non-condensate density 
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via 

^(^^ = / {2^nf e^e>.P) - 1 ' 

as discussed further in [10, 65]. Predictions of the Hartree-Fock theory for condensate 
fraction and atomic density profiles are shown in Fig. 2, showing clearly the separation 
into condensate (dot-dashed, bhic) and thermal cloud (dashed, red). In the absence of 
interactions, the thermal cloud in a harmonic trap acquires a gaussian density profile; 
however, the presence of effectively repulsive atomic interactions between the atoms, 
and thus also between the condensate and the thermal cloud, leads to the appearance 
of a local dip in the thermal density at the centre of the trap, where the condensate 
has its maximum density. 

The Hartrce-Fock basis will be used in our subsequent analysis (Sec. 4.4.3) as the 
suitable basis for developing more advanced dynamical treatments. Before doing so, 
however, let us first present various alternative (static) formulations, which are more 
advanced in the sense that they include the dressing of particles into quasiparticles, 
and thus appropriately generalize the Bogoliubov equations of Sec. 2.2. 



3.1.3. The Hartree-Fock- Bogoliubov (HFB) Limit: This is a generalization of the 
Hartree-Fock regime, in which one includes all quadratic non-condensate operators 
to the hamiltonian, and also explicitly maintains correlations of pairs of like non- 
condensate operators, i.e. the anomalous average r7i(r,t) = (5(r, t)(5(r, t)) and its 
conjugate rh*(v,t) = {5'^ {v,t)5^ {v,t)). Inclusion of these additional contributions 
{5H^°^ + 5H^^° + 5Hi^^) reduces the system hamiltonian to the following 
quadratic form, known as the Hartree-Fock-Bogoliubov hamiltonian: 



H « Hhfb = {Ho + SHo) + (^1 + + (^2 + SH2) . 



(64) 



This effective hamiltonian can be diagonalized via a Bogoliubov transformation of 
Eq. (47). Using the splitting (j) = (pQ S introduced earlier, the grand canonical 
Hamiltonians (45) and (46) take the following generalized 'mean field' forms 

k[= j dr ^(ho + g\(j)o\^ + 2gho - f.^ (1)0 + gmo(j)*Q +h.c., (65) 



S5 (<^; 



) ^^^^] } 



(66) 



As before, the equation describing the condensate is obtained by the requirement that 
K[ = 0, which leads to the following time-independent form of the generalized GPE 



ho + g\(l)of + "^gfio <Po + gmo(lio = MfAo • 



(67) 



(68) 



Diagonalizing the K'2 term (66) with the Bogoliubov transformation (47) leads to the 
generalized Bogoliubov equations (c.f. Eqs. (40)): 

£(r) M{r) \f u,{r) \ / u,{r) 
-M*{r) -£*(r) J \ Vi{r) J [ Vi{r) 

where 

£(r) = L(r) + 2gfio{r) = ho{r) + 2g\M^)\^ + 2gM^) " M (69) 
Mir) = M(r) + gmo{r) = glM^^)]^ + gmoiv) . (70) 
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Thus, the excitations of the system are now BogoUubov quasiparticles, whose energies 
arc obtained from Eq. (68). At finite temperature the Bogoliubov excitations will be 
thermally populated. Since we are considering a static thermal cloud, the occupation 
of the quasiparticle states is diagonal, i.e. 0l(ij)o = ^ijfi and 0iPj)o = 0. Here we 
have defined the so-called quasiparticle distribution function /j, defined by 

where the subscript denotes a static value. Thus, the equilibrium values of normal 
and anomalous averages can be expressed in terms of these distribution functions and 
Bogoliubov amplitudes Ui{r) and Vi(r) via [66] 

no(r) = {6Hr)S{r)) = ^(|«,(r)p + Mr)\')f^ + \vi{r)\', (72) 

i 

mo(r) = {S{r)S{r)} = J2^,{r)v*{r){l + 2/°). (73) 

i 

The set of equations (67)- (73) are collectively termed the time- independent 
Hartree-Fock-Bogoliubov (HFB) equations. They can, in principle, be solved self- 
consistently [69]. Hence the energies and eigenfunctions of the quasiparticles can 
be found , as well as the equilibrium condensate, 'normal' and 'anomalous' averages 
(j^cO = \4'o\'^; "0 and toq respectively) - see e.g [66] and references therein. 

The HFB equations can be easily seen to reduce to the HF theory in the 
appropriate limit. To visualize this, consider setting the quasiparticle operators bi 
to be equal to the single-particle operators Sj, which can be achieved by setting 
Vi{r) = in Eq. (47). From Eq. (73) we then see that mo(r) = 0, which implies 
that Ui{r) reduces to a (dressed) single-particle eigenstate. The Bogoliubov equations 
thus reduce to 

£(r)ui(r) = eiUi(r) (ho + 2g\(l)o\'^ + 2gn<^ Ui{r) = eiUj(r) , (74) 

thus mapping the excitation energies onto the HF energies of Eq. (61). 

From one perspective, the inclusion of the anomalous average within HFB is 
desired, as it provides a better approximation to the (many-body) wavefunction, 
and hence it lowers the free energy of the system|l. However, as formulated above, 
these equations have two fundamental limitations which restrict their applicability. 
Firstly, the anomalous average defined by Eq. (73) diverges as i — > oo, which is 
a direct consequence of our somewhat careless handling of the effective interaction, 
thus requiring us to revisit the imposed pseudopotential approximation of Eq. (16). 
Secondly, the homogeneous spectrum of elementary excitations predicted by the HFB 
equations does not vanish, as it should, in the zero momentum limit, i.e. it exhibits 
a 'gap'. Both these subtle issues can be circumvented by more careful considerations, 
as discussed in the subsequent section. 

3.2. Introduction of an Effective Interaction: 

In our discussion so fax, we have simply replaced the interatomic potential V{r — r') by 
an effective contact potential gS{r — r'), via Eq. (16), where 5 is an effective constant 

II Another appealing feature of the HFB approach is that it is 'conserving', in the sense that the 
response functions generated by this approach can be shown to satisfy various conservation laws. A 
discussion of conserving vs. gapless approximations can be found in the seminal paper by Hohenberg 
and Martin [48] , with more recent summative discussions made by Griffin [67] , and additional related 
insight presented in the reviews of Andersen [68] and Yukalov [70]. 
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strength of binary interactions; as this approximation is actually related to the above- 
mentioned problems arising with the HFB formulation, we should now justify its use 
and discuss how to overcome any associated limitations. 

It is natural to wonder in what sense the exact non-local interatomic potential 
V{v — r') between two atoms can be meaningfully replaced by a contact potential. 
The answer actually lies in the details of the particular system under consideration: 
For the low temperature dilute systems we are considering here^, atoms spend most 
of their time far apart from each other, i.e. at distances much larger than the typical 
range of the interatomic potential, such that short-range correlations are unimportant. 
Typical atomic interactions can thus for most practical purposes be described as 
scattering processes, and one is only interested in the effect that the full potential 
has at large distances |r — r'|, i.e. in the asymptotic scattering states. The only effect 
of atomic interactions on these states is a change of phase of the quantum-mechanical 
wavefunction [10, 31, 71, 72]. 

It turns out that such a phase shift can be well-reproduced by means of a so-called 
'pseudopotential'. To explain this, we consider the scattering of two atoms located at 
a relative distance r from each other, where r is much larger than the effective range of 
their interatomic potential. At sufficiently low energies, only one scattering channel is 
energetically accessible, the so-called s-wave scattering one. In this case, the net effect 
of the potential can be well-approximated by a pseudopotential of the form [31, 35] 

d 

Vpseudo {r) = 96{r) — r g5 (r) . (75) 

It is important to note that the above pseudopotential is actually an operator which 

acts on the particle wavefunction, thus leading to terms of the form Vpsoudo(^)'/'(^) = 
g5{r)d{r(j))/dr. However, the problem simplifies considerably under the assumption 
that the pseudopotential acts on unperturbed free-particle wavef unctions; in this limit, 
the operator (d/dr) can simply be replaced by 1, reducing the problem as indicated 
by the arrow in Eq. (75). However, extreme care is needed when using a hamiltonian 
based on Eq. (75), and there are in particular two related issues which should be 
considered: firstly, under what conditions can the presence of the operator d/dr be 
overlooked, and secondly what is a suitable value (or function) for g which would then 
correctly describes the scattering process? 

To address both of these, we note that a scattering process can be viewed as the 
end result of a collisional process which in general includes the effect of repeated virtual 
collisions between the two atoms; mathematically, this corresponds to the repeated 
action of the exact interatomic potential V{r). While the effect of each repeated 
collision may be very laxge, one can actually construct an infinite series over such 
virtual collisions, i.e. a series in powers of the exact interatomic potential and their 
corresponding propagators. This series actually converges, and such a mathematical 
construction amounts physically to upgrading the exact interatomic potential to an 
effective one, known as the Transition (or simply T) matrix; the defining relation of 
the T-matrix in terms of the actual interatomic potential, the so-called Lippmann- 
Schwinger equation [10, 73], is depicted diagrammatically in Fig. 1 (bottom). It is 
important to realize that the pseudopotential approximation gS{r) is only strictly 
meaningful when imposed onto this upgraded effective interaction. 

^ Typica l conditions fo r this are na^ -C 1 and a -C XdB, where a is the s-wave scattering length and 
XdB = \/2nh^ /mkBT is the thermal de Broglie wavelength of the atoms. 
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So, what would happen if one were to overlook this, and replace the exact 
interatomic potential in the full system hamiltonian of Eq. (15) by g5{v — r') as 
hinted by Eq. (16)? As long as the subsequent treatment is restricted to first 
order perturbation theory, this would be fine. However, imagine one chose instead 
to use this hamiltonian to compute properties of the system within second order 
perturbation theory. Considering for simplicity the case of a homogeneous Bose gas, 
second order perturbation theory would generate terms containing the interaction 
strength contribution g*^^) = g[l + fj^pi^/Sp)] where the prime implies that the 
summation is restricted to nonzero momenta - see [74] for more details. However, 
this contribution is easily seen to diverge for large momenta p (so-called ultraviolet 
divergence). Such a divergence has appeared here because, without realizing it, 
the above procedure has resulted in double-counting (due to the use of a constant 
interaction strength combined with an unrestricted summation over momenta). The 
reason is that, in this second order perturbation theory, the effective interaction 
strength g should appear for the total interaction strength to second order, and 
not for each of its individual contributions as implicitly assumed in the above g^"^^ 
expression. In this context (i.e. when the pseudopotential approximation has already 
been made in the original hamiltonian), the problem can nonetheless be solved by the 
technique of 'renormalization' whereby one simply replaces the second-order expression 
for the effective interaction strength g^^^ by g [75, 74]. However, one can avoid the 
need for such renormalization by ensuring that the pseudopotential approximation 
is correctly implemented in the calculations; in other words, the pseudopotential 
approximation should only be made after the exact interatomic potential has been 
upgraded to the effective interaction given by the T-matrix. In fact the problem 
above arose because the pseudopotential approximation was imposed on the first order 
contribution (Born approximation) in the series expansion of the T-matrix in terms 
of the exact interatomic potential. 

We should still address the related issue of the strength of the effective interaction 
given by the T-matrix; in general, this should depend on the momentum of the 
incoming particles [73, 76]. For very dilute g natural first assumption would 

be that the scattering process is taking place in vacuum; in that case, the resulting 
effective interaction between two atoms is known as the two-body T-matrix. As long 
as one is only considering atoms whose relative momentum is small, one can ignore 
the momentum dependence of the T-matrix, thus approximating it by a constant, g, 
equal to the zero-energy, zero-momentum limit of the full T-matrix. The value of this 
latter quantity is fixed by a single parameter, namely the shift of the projection of the 
asymptotic wavefunction onto the relative coordinate axis from the origin - a quantity 
known as the s-wave scattering length, a. The scattering length can be either positive, 
or negative, depending on the intricate details of the exact interatomic potential, with 
the respective sign indicating effective repulsion, or attraction, between the atoms. 
It is important to note that this constant value is only a good approximation when 
describing collisions between low-momentum states, and it can be thought to arise by 
the elimination of high-lying modes from the problem [77, 78]. In fact, for momenta 
larger than h/a, the full two-body T-matrix would rapidly decrease to zero, an effect 
not reflected in our choice of a constant strength g, which would thus need to be 
implemented in conjunction with an upper momentum cut-off [36, 37, 49]. 

We have thus discussed how a contact interaction potential arises in the context 
of a microscopic theory of ultracold Bose gases, which is actually a crucial point in 
assessing the validity of any ab initio theory. As long as this replacement can be 
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routinely performed, one does not need to worry about the above subtle issues, i.e. 
one can start off from the hamiltonian of Eq. (17); this is, for example, true in arriving 
at the HF theory of Sec. 3.1.2. However, in trying to construct more advanced theories 
- such as the HFB - from a microscopic perspective, one should be more cautious, and 
investigate the following issues, as expounded clearly by Burnett [79]: 

(i) Is the introduction of the rekn'ant effective interaction (T- matrix), which is then 
replaced by a contact pseudopotential, performed consistently at the level of the 
considered approximation? 

(ii) How does the fact that collisions in a BEC take place within a medium, rather 
than in vacuum, affect the above considerations? 

3.3. Problems of the HFB Formulation: 

HFB differs from the simpler HF by the explicit consideration of the static anomalous 

average, TOo(r) of Eq. (73). Since the HFB equations, Eqs. (67)-(73), are formulated 
in terms of the constant interaction strength g, the summation in Eq. (73) should 
exclude high-lying modes, as these have already been implicitly considered in replacing 
V{v — y') by the effective interaction in vacuum (two-body T-matrix). If this were not 
done, then the resulting expression for 77io(r) would blow-up as i ^ oo, i.e. it would 
exhibit an ultraviolet divergence. Having identified the origin of this divergence, we can 
routinely eliminate it from the problem, by regularizing the anomalous average. This 
is achieved via the subtraction of the high energy part, whose contribution corresponds 
precisely to the difference between the actual and the effective interatomic potential 
(which would otherwise be double-counted). We thus write [79, 78] 

mo(r) m^{r) = mo(r) - lim,^^Ui{r)v* {r) . (76) 

Any subsequent treatment including the anomalous average has to be performed 
with this regularized form, so we will henceforth always quote mQ^(r) as the 
anomalous average, implicitly assuming that such an essential 'renormalization' has 
been performed. 

The second related problem arising from the inclusion of the anomalous average 
is the appearance of a gap in the excitation spectrum of the homogeneous gas. 
As already done for T = in Sec. 2.2, we consider here the predictions of the 
HFB equations in the case of a uniform gas, for which T4xt(r) = and 0o(r), 
no(r) and m^(r) are constant. From Eq. (67) the chemical potential becomes 
/u = g(|0oP + 2no + itiq). On the other hand, the energy spectrum is given by 
e{p) = [p^ /2m + 2g(|(^oP + ^o) ~ /^]^ - .9^(100^ + m^)^. The dispersion relation 
for the quasiparticles is obtained by substituting the value of /x into this expression. 
However, for p ^ the energy spectrum has a finite value (proportional to mg*), or in 
other words, it has a gap. This is however in direct contradiction with the Goldstone 
theorem; the latter requires that the energy spectrum arising from a theory based 
on symmetry-breaking should be gapless, which is equivalent to the statement that it 
should cost zero energy to excite the lowest mode (known as the Goldstone mode) [80] . 
The HFB theory must thus be modified to be in accordance with this theorem. The 
simplest, yet somewhat heuristic, way to solve this problem is to completely ignore 
the 'anomalous average' ifiQ. Such an approximation was routinely made in the early 
literature, and is often called the HFB-Popov approximation, as discussed by Griffin 
in [67, 81], although objections have been noted as to the use of this terminology [82]. 
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3.3.1. The HFB-Popov (HFBP) Limit: The Hartree-Fock-Bogoliubov-Popov 

(HFBP) limit is an intermediate regime between Hartree-Fock and Hartree-Fock- 
Bogoliubov. In comparison to HF, it additionally includes the contribution ff^"^*^ 
which dresses the finite-temperature single-particle energies to quasiparticle energies, 
whereas it can also be obtained from the full HFB theory by discarding all 5Hf~"^ 
contributions [i = 0, 1.2). Mathematically, this can be expressed as 

TJ £71 ttBOG 

^HFBP — J^HF + J^2 

= Hhfb {SHr'' + SHr"" + SHi'''') , (77) 

with the system liamiltonian in this regime given by 

Hhfbp = {Ho + SH^^) + + dH^) + + 5H^^) . (78) 

This leads to the same GPE as in the HF limit, i.e. the condensate is described by Eq. 
(60). However, contrary to the HF limit, quasiparticle dressing is actually included 
here via the Bogohubov equations of Eq. (68), in which the respective operators are 
now modified according to: 

£(r) ^ £p(r) = £(r) = L(r) + 25no(r) 

= /io(r) + 2.g|(/)o(r)P + 2.gno(r) - H (79) 
M{v) ^ Mp{y) = M{v) - gmo{v) = M(r) = g[MT^)f . (80) 

It is easy to verify that this set of equations leads to a gapless energy spectrum 
(and hence corresponds to a better theory for calculating frequencies of elementary 
excitations). However, this approximation may also be problematic at very low 
temperatures T <C T^, since then mo is of the same order as no [82]. 

Takano [41] was, to the best of our knowledge, the first person to point out 
(using slightly different arguments to those presented above) that the gap problem can 
actually be removed by going beyond the mean- field approximations of Eqs. (51), (53) 
via the consistent consideration of correlations of three fluctuation operators. The 
next section presents two possible (yet somewhat heuristic) generalized approaches 
along those lines. 

3.4- Static Generalized Many-Body Theories: 

Having identified the difficulties with the HFB model, one can now attempt to 
construct an improved theory which includes the effect of the pair anomalous average, 
but simultaneously also yields a gapless homogeneous excitation spectrum. 

The strength of the effective interaction employed thus far corresponds to 
collisions taking place in vacuum, whereas in our system there is actually an active 
atomic medium present. This will in general affect the scattering in two ways [79]: (i) 
firstly, the intermediate states - denoted by thick (blue) lines and filled arrows in Fig. 
1 - may actually be thermally occupied, thus providing a bosonic enhancement for the 
transfer rates into those states; (ii) moreover, the intermediate states may actually be 
dressed (quasiparticle) states, instead of single-particle ones. The anomalous average 
mo(r) contains information about correlations between two nearby atoms. We have 
already argued that its contribution over high-lying modes is related to the difference 
between the actual and the effective interatomic potential. It should therefore come as 
no surprise that its static value, ih^, over low-lying modes contains information about 
modifications to the effective interaction due to the presence of a medium; in technical 
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jargon, inclusion of ?no(r) has the effect of upgrading the previously-considered two- 
body T-matrix describing scattering in vacuum to the so-called many-body T-matrix 
[73, 81], which describes the additional modifications due to the presence of the 
medium (see also [83, 84] for a more detailed handling of many-body theories which 
also includes the regimes of lower dimensionality); these theories thus clearly extend 
beyond the HFB limit. 

Following Proukakis, Morgan, Choi and Burnett, one can thus define a generalized 
effective interaction as [85]: 

pW=.(l+f). (81) 

Such a definition enables us to express the interaction term between two condensate 
atoms in the form g(r)|c6op0o = .9 [I'/'oP + to^] (/io- A careful consideration of the 
HFB equations reveals that such a generalized interaction does not enter in all terms 
of the HFB equations; in particular, even though g{v) appears explicitly in the finite 
temperature GPE of Eq. (67) , there is nonetheless no ih^ contribution within the C. 
operator of Eq. (69). This leads to an inconsistency in the manner in which interactions 
are handled within the HFB formalism, and can be identified as the cause for the 
appearance of a gap in the homogeneous excitation spectrum at low momenta. This 
problem is a direct consequence of the truncation of the coupled equations of motion 
to some particular order [68], and has already been identified in [86]. Nonetheless, 
it has been suggested [85] that this problem can be 'fixed by hand', by forcing the 
interaction strength in the HFB equations to take the modified form given in Eq. (81). 

Although the above arguments hold for the collisions between two condensate 
atoms, an immediate issue arises regarding the form of the interaction strength when 
one of the atoms belongs to the thermal cloud. Consistency requires this interaction 
to be also computed in the many-body limit. However, due to the large range of values 
of the relative momenta between the two colliding atoms, a full treatment, expounded 
analytically in [73, 76], is numerically quite involved, and one typically resorts to 
approximations. In fact, an interpretation of the generalized effective interaction 
g{v) can be made by resorting to the homogeneous limit where the role of many- 
body effects is well-known. Such a comparison reveals that the generalized effective 
interaction of Eq. (81) indeed includes one aspect of the many-body effects, namely the 
population of intermediate states, but it does not include the dressing of these states 
to quasiparticles. Formally, this is equivalent to the zero-energy zero-momentum limit 
of the full (many-body) T-matrbc. In this case, the many-body T-matrbc describing 
the collision between a condensate and a thermal atom reduces to the corresponding 
two-body expression in vacuum, i.e. the interaction strength is simply equal to g. In 
the other extreme, one can approximate the effective interaction strength between a 
condensate and a thermal atom by the same expression as for two condensate atoms, 
i.e. £f(r) - the latter expression was also subsequently derived by rigorous consideration 
of pseudopotentials by Olshanii and Pricoupenko [87]. 

The above physically-motivated, yet mathematically ad hoc procedure leads in 
general to the following somewhat generalized set of Bogoliubov equations 

/ C'{v) M'{v) \ ( u,{v) \_( u,{v) \ 
\ -M'\r) -C'\r) ) [ v,{r) J ' { v,{r) J ' 

where 

£'(r) = ho{r) + 25e(r)|0o(r)|' + 25*(r)no(r) - /x (83) 
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Position in trap Pooilion in trap 

Figure 3. Equilibrium spatial dependence (near the trap centre) of the 
generalized mean fields for a one-dimensional trapped gas consisting of 2000 
atoms at a temperature T/Tc ~ 0.6, where there is an approximate condensate 
fraction of 50% (here denotes an approximate ideal gas degeneracy temperature 
[96]); densities are plotted in units of l^^ and position in units of Iz- Top 
Row: Condensate |</>o(2;)P (left) and thermal cloud no(2) (right) computed self- 
consistently for the generalized HFB approaches (pc = 3(2)) with gt = g(z) 
(dashed) or gt = g (solid) and for HFB-Popov (dot-dashed); dotted lines 
indicate corresponding results in the limit when interactions between condensate 
atoms are ignored, and only interactions between condensate and thermal atoms 
are maintained. Bottom Row: Corresponding profile for the pair anomalous 
average mQ-{z) within the two generalized HFB approaches (left); Typical spatial 
dependence of the generalized effective interaction g{z) within the two generalized 
HFB theories, plotted in units of the constant strength g, at different temperatures 
- from top to bottom: T/Tc 0.1, 0.2, 0.6 and 0.9. (Reprinted figures with 
permission from N.P. Proukakis, S.A. Morgan, S. Choi and K. Burnett, Phys. 
Rev. A 58, 2435 (1998). Copyright (1998) by the American Physical Society.) 



M'{r) = QelMr)? 



(84) 



Here .9c(r) denotes the interaction strength between two condensate atoms given by 
9c{^) = ^(r), whereas 3t(r) expresses the effective interaction strength between a 
condensate and a thermal atom. Clearly these equations are also coupled to a suitably- 
generalized GPE, given by 



ho + fi(c(r)|0o|^ + 25ft (r) no (r) (^o(r) = M0o(r) 



(85) 



The above discussion leads naturally to two different theories, which we shall 
henceforth collectively refer to as generalized HFB, depending on whether (/((r) = 
£(c(r) = g{r), or S(t(r) = g . Clearly the limiting case gd'r) = fl't(i') = g corresponds 
preciselt to the HFB-Popov limit discussed in Sec. 3.3.1. Such generalized theories were 
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applied with some success to the study of finite temperature excitation frequencies [90] 
(see Sec. 7.1), to vortices [89], and to the coherence of 2D condensates [91, 92]. Note 
that a related approach for the self-consistent calculation of the coupling constant 
appearing in the GPE was suggested in [93], whereas alternative approaches to 
overcome the problem of the gap of the HFB theory arc discussed in [82, 94, 95]. 

It is appropriate to make a brief comparison of the approaches mentioned above, 
as discussed in [85] (see also [90, 66, 89]). Fig. 3 displays the position dependence of 
the generalized mean fields for the previously discussed theories. Focusing initially on 
the normal averages 0o and fio, we note that the predictions of all such theories differ 
considerably from the corresponding ones in the limit when condensate-condensate 
interactions arc completely ignored, and the only interactions considered arc those 
between condensate and thermal atoms (in the various approximations); in particular, 
inclusion of repulsive condensate-condensate interactions leads to the condensate 
wavcfunction being more spread out and having a lower central density; in comparison, 
the differences in the respective profiles predicted by the various interacting theories 
are generally small, with the only noticeable deviation arising for the case gt = fl'(r). 
Investigation of the change in the temperature dependence of the condensate fraction 
again reveals marginal differences, with the corresponding predictions for excitation 
frequencies only being significantly different for the case = g{r) (see also Sec. 7.1). 
The anomalous average fh^ is found to be negative and to exhibit minima at opposite 
off-centred points, with its value exhibiting a local maximum at the trap centre. As a 
result, the generalized effective interaction ^(r) is consistently smaller than, or equal 
to, g. In particular, g(r) exhibits a local maximum at the centre of the trap (where 
the condensate density is largest) and approaches the value of g asymptotically far 
from the trap centre; moreover, its central value g{Q) is temperature dependent and 
reaches its minimum value at the transition temperature. 

3.5. Brief Summary 

The mean field effect of a static thermal cloud on the condensate can be calculated by 
an approximate inclusion of the {H^ + H^) contributions to the hamiltonian by means 
of the mean field approximations of Eqs. (51) and (53). These approximations give rise 
to a set of generalized mean fields consisting of the condensate mean field, <j){v), the 
normal average denoting the non- condensate or thermal density, n(r) = (5'l'(r)^(r)), 
and the anomalous average, Tfi{v) = {S(r)S(r)). Implementation of Eqs. (51) and 
(53) thus reduces the system hamiltonian to an approximate quadratic form (more 
general than that considered in Sec. 2) which enables the self- consistent evaluation of 
(the static values of) n(r) and m(r) in term,s of the Bogoliubov functions Ui(r) and 
Vi{r) (via Eqs. (72)- (76)). These approximations cannot describe dynamical effects as 
they discard contributions describing both particle-exchange collisions between thermal 
and condensate atoms and collisions between pairs of thermal atoms which lead to 
thermal population redistribution. There are four variants of such generalized mean 
field theories, as summarized below and in Table 1. Note that an alternative review of 
(most of) these theories using the Green's function formalism, can he found in [81]. 

• Hartree-Fock (Eqs. (62)-(63)): Only maintains quadratic hamiltonian contribu- 
tions of the form. Hhf = J dr6^{---)6, i.e. contributions containing one creation 
and one annihilation non-condensate operator. Usually implemented in the semi- 
classical approximation. 
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• Hartree-Fock-Bogoliubov (Eqs. (67)-(70)): Additionally includes terms of the 
form J dr{- ■ ■)S6 + h.c, i.e. single-particle states are dressed to Bogoliuhov 
quasiparticles. Produces lowest ground state energy, but suffers from an 
inconsistent treatment of atomic collisions, associated with the incorrect 
introduction of an effective interaction which does not treat all collisions 
consistently. 

Two modified mean-field approaches have been proposed to address this issue and 
generate consistent theories in different levels of approximation: 

• Hartree-Fock-Bogoliubov-Popov (Eqs. (60), (68), (79)-(80)): Contributions of the 
form J dr{- ■ •)M + h.c. are only maintained in the limit mo{r) = - often referred 

to as the 'Popov ' approximation. 

• Generalized Hartree-Fock-Bogoliubov (Eqs. (82)-(85)): The T > Bogoliubov 
equations are amended 'by hand' to include physics beyond the usual HFB 
quadratic hamiltonian, such that the effective collisional strength between two 
atoms becomes upgraded throughout the coupled equations from g to g{r) = 
g{l + mQ-/0§) for condensate-condensate collisions; a variant of this theory also 
includes a similar replacement for collisions involving one condensate and one 
thermal atom. 

Having concluded our discussion of static theories, we next extend our formalism 
to incorporate effects introduced by the dynamics of the thermal cloud. 



4. Finite Temperature Mean Field Theory: Dynamic Case 

Our discussion so far has been restricted to the study of static variables, aimed 
mainly at interpreting finite temperature properties of the ultracold gas, such as 
density profiles and coherence properties at steady state. This treatment can be 
systematically generalized to time-dependent variables. We start by obtaining exact 
coupled equations of motion for the condensate and the non-condensate components, 
and then (lisc;uss various approximate treatments which enable direct numerical 
simulation of the full coupled evolution. 



4-1- Exact Dynamical Evolution 

Using the Heisenberg equation of motion, Eq. (19) we readily obtain the following 
exact evolution for the condensate mean field [97, 88]: 

-dHr,t) 



= {ho^{r, t)) + (r, t)^{r, t)^{r, t)) . (86) 

Using the symmetry-breaking decomposition into condensate and non-condensate 
contributions via Eq. (22), and suppressing explicit spatial dependence for brevity, 
the latter term can be re-expressed as 



(87) 



This then yields 



ho + g\4>\^ (f + + + g{SUS) 
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Clearly the condensate evolution is dynamically coupled to the behaviour of the 
non- condensate via the evolution of normal and anomalous averages. Using similar 
arguments, we obtain the exact expression 



ih 



d 5{v, t) 

dt 



5(v,t),H]=ih^^ (*(r,t)-(*(r,i))) 



(89) 



which can be written as [109] 



ih 



dS 
di 



hoS + g 
+ g (S'^U - 



+ 



-{5^5))+g(l>* (55 -{U)) 



(90) 



One can thus obtain the exact equations for the evolution of normal and anomalous 
averages either directly from Eq. (90), or equivalently via 

d . 



66, H 



(91) 



While formally exact, without any truncation this procedure leads to an infinite 
hierarchy of coupled equations of motion for the higher order correlations. This is easy 
to see, since for the given hamiltonian, the right-hand side of the equation of motion 
for a product of n operators actually contains n + 2 operators. Thus, in order to solve 
such equations one must choose an appropriate set of generalized mean fields which are 
assumed to accurately describe the system to the desired level of approximation. This 
Section discusses the various approximations that are typically made, and the physics 
that they contain. An alternative, yet closely related approach based on the method 
of non-commutative cumulants and pursued by Kohler, Burnett and collaborators is 
discussed in Appendix B. 

In order to formulate such theories, we occasionally need to resort to a basis 
set notation in terms of single-particle eigenstates ipi{r). Following the notation 
introduced in the book of Blaizot and Ripka [69] and implemented extensively in the 
group of Keith Burnett (Oxford) the condensate and non-condensate contributions 
are denoted by 



i 



Zi = (ai) , 



{a-i) . 



Correspondingly, the normal and anomalous averages are defined by 



Pij = {cUi) 



{CjCi) . 



(92) 



(93) 



For convenience in our subsequent discussion we also define the 'Hartree-Fock 
hamiltonian' (h) and the 'Pairing Field' (A), along with their 'reduced forms' h" 
and A"^, via their respective matrix elements 



K = h-j + J2 VikijZ*kZi = {i\ho\j) + 2 ^ V.kji {zlzi + pik) , (94) 

kl kl 
^ij = ^Ij + X! ^ijklZkZl = ^ Vijkl (ZkZl + Klk) . (95) 



kl 



kl 
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Here, Vikji = {ik\V\jl) = {l/2)[{ik\V\jl) + {ik\V\lj)] is the symmetrized form of the 
interaction matrix clement defined by Eq. (11), corresponding to the collision of two 
atoms initially in states j and I, which emerge from the collision in states i and k; we 
also note here that h*j = hji. 



4-. 2. The Hartree-Fock Approximation 

As in Sec. 3.1.2, the simplest finite temperature approximation that one can consider 
is formulated in terms of the condensate ^(r, t) and the thermal n(r, t) components. In 
the position representation, Eq. (88) leads to a generalized finite temperature Gross- 
Pitaevskii equation 

'^^^di^ = [ho + g{\Hr,t)\^ + Mr,t))] cl>{r,t) , (96) 

which is the time-dependent generalization of Eq. (60). 

In order to express the coupled equations for condensate and non-condensate, we 
shift to our basis set notation, where Eqs. (88) and (90) generate the following set of 
self-consistent equations [69]: 

ih^ = h'^z + A^z* , ih^ = [h,p] , (97) 

which constitute the time-dependent generalization of the 'Hartree-Fock' equations of 
Sec. 3.1.2. As before, this level of approximation only inc-ludes mean field coupling 
between the two subsystems and will not be considered further here. We will however 
return to an appropriate generalization of these equations which additionally includes 
particle exchange in Sec. 4.4.3. 



4-.3. The Hartree-Fock-Bogoliuhov Approximation 

In the Hartree-Fock-Bogoliubov (HFB) approximation, one works with slightly 
generalized mean fields: in addition to the condensate mean field and the normal 
average, one further includes the pair anomalous average rh{v,t), i.e. one assumes the 
system is well-described by the HFB hamiltonian of Eq. (64). 

In first instance, we consider weak perturbations around the equilibrium solutions, 
which is sufficient to yield expressions for the dominant damping mechanisms both at 
zero and finite temperatures. We then proceed to give the full dynamical equations. 



4-3.1. A Perturbative Linear Response Treatment: One assumes that the system is 
close to steady-state, with both the condensate and the thermal cloud described by 
small deviations around their equilibrium values, via 

^(r,t)~e-'''*/'^[0o(r)+^0(r,t)] , (98) 

n(r, t) ~ no(r) + dn{r, t) . (99) 

This is a generalization of the linear approximation used to derive the (zero- 
temperature) Bogoliubov equations in Sec. 2.2. Consistency requires one to also look 
at small variations about the anomalous average rh{r,t), via 

?fi(r,t)~e-2*^*/''[m^(r) + (5m(r,t)] . (100) 

We will see that this approach leads to the inclusion of the dominant damping 
mechanisms for elementary excitations. The wavefunction 0o(r) satisfies the static 
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generalized GPE of Eq. (67). Linearization yields the following equation for the 
condensate perturbations [88, 98, 99, 100] 

ih—5(j){r,t) = [/io-/x + 2fif(|(/.o(r)|2 + flo(r))J 5(j>{r,t) 

+ [g\M^)\\T)+gm^{T)\5cl>*{T,t) 

+ 2#o(r)<5n(r, t) + #o(r)<5m(r, t) . (101) 

As discussed earlier, the static values of no(r) and ■mQ{r) arc obtained via a 
transformation to a quasi-particle basis - see Eqs. (72)- (76). Clearly the condensate 
perturbations are coupled to the dynamic change in the thermal cloud 6h{r), and the 
anomalous average drh{r) which respectively acquire the forms [98, 99, 100] 

6h{r,t) = J2 { K (r)«,-(r) + vUr)vj{r)] f,j{t) 

ij 

6rh{r,t) = ^ |2<(r)uj(r)/ij(t) + Ui{r)uj{r)g,j(t) 

ij 

+ v:irX^ir)g;.{t)}. (103) 

Here both fij = {filPj){t) — Sijff and gij correspond to non-equilibrium distribution 
functions (/° denotes the Bose-Einstein distribution of Eq. (71)). The evolution of 
such functions can be easily obtained within our generalized quadratic (quasiparticle) 
hamiltonian. Here wc apply perturbation theory to lowest order, i.e. integrate the 
equation of motion for fij and gij which yields corresponding expressions to lowest 
order in the interaction strength g. Any contributions of Sh(t), and 6rh{t) appearing 
on their respective right-hand sides, will thus, to this order, yield no contribution to the 
expressions of interest. With that in mind, the evolution of the above non-equilibrium 
distribution functions is governed by the simplified equations: 

d_ 
"dt- 



i^irJij it) = i^j - ^i)fij (*) + ^gAij (/° - /°) , (104) 



ifi§^9ij it) = {cj + ei)fij it) + 2gBij (1 + + /?) , (105) 

where we have dropped from their respective right hand sides all terms containing 5n 
and 5fh. The coefficients Aij and Bij appearing here are sums of overlap integrals of 
the condensate with the respective quasi-particle amplitudes, given by [99, 100] 

+ 5(j>* {uiU* + ViV* + UiV*) } (106) 

+ sr («•»• + <»• + <.'«;) } . (107) 

Importantly we find that the above procedure already incorporates the two 
main damping mechanisms encountered in experiments with ultracold gases. To 
visualize this, we consider oscillations of the condensate at some given frequency w. 
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i.e. 5(f){v, t) = (5</)o(r)e~*"*. These condensate fluctuations induce changes in fij{t) and 
gij{t), which modify Sh{r, t) and Srh{r, t) and in turn act back onto the condensate via 
Eq. (101). This leads both to a shift of the frequency of oscillation and to damping; 
the damping is determined by the coefficient 7 = 71, +7b, which consists of two terms: 

(i) A thermal excitation (quasiparticle) of energy can interact with the collective 
mode of the condensate of typical energy fiwo, and thereby become converted into a 
higher energy excitation, Cj — ei + fujo (i.e. energy is absorbed from the condensate). 
Such a process (Landau damping) damps the condensate oscillations at a rate 

7L = W E l^^i I' ifi - fi) ^ (^0 + ' (108) 

where the factor {f° - /«) = [/°(/° + 1) - (/« + 1)/°] denotes the difference in the 
amplitudes for the destruction of a quasi-particle of energy Cj, and the simultaneous 
creation of a quasi-particle of energy ej and its inverse process. 

(ii) The absorption of a quantum of oscillation can also lead to the creation of two 
excitations, of respective energies Cj and ej with gj + Cj = fkJo- This process (Beliaev 
damping) leads to damping at a rate 

7B = 2V E l^'^l' (1 + /° + fj) ^ (^0 -ei- ej) , (109) 

where the factor (1 + ff + /«) = [(/« + l)(/j' + 1) - /O/o] denotes the difference in 
the amplitudes for the simultaneous creation of two quasi-particles of energies e, and 
Cj and its inverse process. 

In the homogeneous limit, the 'low-temperature' (quantum) regime fcsT <^ fUvo 
is dominated by Beliaev damping, while the 'high-temperature' regime ksT ^ hu)o is 
dominated by Landau damping [78, 99, 100, 101, 102]. 

4-3.2. Self- Consistent Time-dependent Hartree-Fock-Bogoliubov Theory: Rather 
than making a linear response analysis for situations close to equilibrium, one can 
formulate the Hartree-Fock-Bogoliubov problem in a fully coupled dynamical manner. 

We are thus interested in writing down full dynamical equations for the selected 
mean field variables cf){r,t) = {cf){r,t)), h{r) = {S^{r,t)S{r,t)), and 7fi^{r) = 
((5(r, t)6{r, t)). These can be obtained from Eqs. (86), (89) and (91) by replacing H by 
Hhfb- Alternatively, one arrives at the same equations using the full expressions of 
Eqs. (88) and (90) and subsequently imposing the decoupling approximations of Eqs. 
(51)-(53), in order to reduce the expressions to a closed system of equations. Note 
that in the HFB approximation {5^66) = 0, so that the equation for the condensate 
evolution differs from the exact equation, Eq. (88), by the absence of the last term. 

There are various ways to write down the precise form of the HFB equations. It is 
convenient to generalize the notation of Eqs. (92)- (95) by defining generalized density 
matrices for the condensed (Re) and uncondensed {Rnc) parts of the system by 

= ( ) ' ^'^'^ ^ ( (P* + 1) ) 

where 1 is the unit matrix. The corresponding generalized hamiltonians are given by 

Hc=(^_ (^c)* _ ^c)* ) , H^c=(^_ -th)*)- ^^^^^ 

thus casting the equations in a notation which is closely related to the conventional 
Green's functions one, used for example in [76, 81]. 
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The time-dependent HFB equations thus take the compact form [111] 

= HcRc and = HncRnc - RncH'^nc ■ (H^) 

The above expressions are a generaUzation of the Hartree-Fock equations of Sec. 4.2, 
and additionally include the damping processes discussed in Sec. 4.3.1. Although 
they contain the essential quasiparticlc physics, these equations have nonetheless been 
formulated in terms of single-particle operators. As a result, their interpretation is 
perhaps not very transparent and one may wish to explicitly perform a Bogoliubov 
transformation to rc-cxprcss these in terms of quasiparticlcs. Such an approach 
was undertaken by Imamovic-Tomasovic and Griffin [103] using an equivalent but 
notationally distinct approach based on the Kadanoff-Baym formalism [43] in terms of 
non-equilibrium Green's functions. This approach is well-documented in the literature 
[43] and we encourage readers who are familiar with the Green's function approach 
to consult Refs. [81] and [104], as well as their extended analysis [103]. We do not 
discuss the HFB equations here any further, as they correspond to a purely mean 
field theory, which does not include the crucial process of particle exchange between 
the condensate and the thermal cloud. This latter issue is addressed in the following 
sections. 

4-4- Dynamical Perturbative Treatment of (^H^ + H^j 

The theories developed thus far include mean-field coupling between the selected 
generalized mean fields. As such they can describe collisions between two condensate 
atoms, or collisions between one condensate and one thermal atom for which the 
number of condensed and non-condensed atoms in the final state remains identical to 
that in the initial state. However, collisional processes which transfer atoms between 
the condensate and the thermal cloud have so far been ignored. In order to include 
these into our treatment, we must thus relax some of the approximations used in our 
previous discussion, essentially following along the lines of important early work by 
Kirkpatrick and Dorfman [105, 106, 107] and Eckern [108]. 

To identify how to proceed, we note that comparison of the exact evolution of 
the condensate mean field of Eq. (88), to its corresponding evolution in the HFB 
basis, reveals the critical absence of the triplet contribution (S^55). This term, whose 
importance has been stressed by one of us (NPP) [97, 88], can in fact be identified as 
describing particle transfer between the condensate and the thermal cloud, a process 
which is however prevented by the mean field approximation of Eq. (53) which dictates 
that (5^56) = 0. Similarly one finds that collisions between two thermal atoms, 
corresponding to the usual scattering processes in a classical gas, can only be described 
if one considers corrections beyond the other mean- field approximation of Eq. (51). 

Collisional dynamics of this form can actually be introduced into the theory 
by perturbatively including the difference between the full system hamiltonian, H, 
of Eq. (23), and an appropriately chosen generalized mean-field hamiltonian -^mf 
[109, 78, 110, 111], i.e. by separating the hamiltonian as 

H = Hmf+(h- Hmf) . (113) 

The first contribution Hmf should be included self-consistently and defines the 
unperturbed basis of the system. This is fixed by the choice of mean field 
approximations (like Eqs. (51), (53)) imposed in reducing {Hs + H4) to approximate 
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quadratic form. For example, setting ^^mf = Hhf defines the unperturbed system 
basis as the Hartrcc-Fock basis, whereas Hmf ~ Hhfb additionally includes 
the pair anomafous average in the basis. Diagonalization of the corresponding 
hamiltonian would define the excitation energies, respectively corresponding in the 
above examples to the dressed single-particle Hartree-Fock, or the Hartree-Fock- 
Bogoliubov quasiparticle energies. 

Such perturbative treatments are presented below for various appropriately 
chosen basis; as the discussion is restricted to second order in the effective interaction 
strength we should comment on the validity of such a methodology: To appreciate 
this, it is important to note that the equations presented here already feature 
a consistent introduction of an effective (T-matrix) interaction; as a result, this 
perturbation theory can effectively be seen as an expansion in terms of the difference 
between scattering in the presence of a medium (namely a partially condensed gas) 
and in vacuum [78]. Thus, although the usual perturbative expansion in terms of 
the actual interatomic interaction would fail, perturbation theory to second order in 
terms of the effective interaction should be valid for a dilute gas sufficiently far from 
the critical region, i.e. when the condition (ksT / gnjVna^ ^ 1 is satisfied [49]; this is 
generally true for experimentally relevant parameters, except in a very narrow window 
near T^. 

The corresponding number-conserving perturbative expansion is presented in Sec. 
5.2, while a related treatment based on the method of non-commutative cumulants, 
which provides a well-defined decorrelation scheme that leads to self-consistent 
expressions in different limits and can also handle 'memory effects' is discussed in 
Appendix B. 



4-4-1- Perturbative Formulation Beyond the HFB Hamiltonian: We start by 
extending the HFB theory of Sec. 4.3.2. CoUisional dynamics can be introduced 
into the theory by the previous perturbative prescription. Using the mean field 
approximations of Eqs. (51), (53) to simplify {H^, + H4) according to Eqs. (54)-(57) 
identifies the mean field 'basis hamiltonian', Hmf as the HFB hamiltonian of Eq. (64). 
In our treatment we thus still focus on the same generalized mean fields 4>, h and m 
as in HFB, but perform second order perturbation theory beyond the HFB basis via 

H = Hhfb + H' . (114) 
where the perturbing hamiltonian, H' , is given by 

H' = H'^ + Hi = [^3 - -^^i] + [-^4 - {SHo + 5H^'\ . (115) 

We thus obtain the exact relations [111] 

in^ = HcRc + Ic, 
at 



cIRnc 
in- 



^HncRnc — RncH\iq^ + Inc ■ (116) 

These expressions differ from their HFB couterparts (Eq. (112)) by the inclusion of 
the 'kinetic' matrices Ic and Inc- These matrices are expressed in terms of averages 
of three and four single-particle fluctuations operators (c^c^^'c) and (c'cP^'cc) beyond 
their approximate mean field values, and are responsible for the inclusion of all relevant 
coUisional terms (detailed expressions can be found in [111]). In principle, one can 
actually obtain a closed system of coupled equations for the evolutions of the matrices 
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Rc, Rnc, Ic and Inc (or equivalently for their corresponding generalized mean 
fields) subject only to appropriate generalized decoupling approximations [97, 111]. 
A potential advantage of such a system of equations would be that it can be solved 
self-consistently using the exact non-local interatomic potential, and does not actually 
require one to resort to the pseudopotential approximation. However, computation 
with actual eigenstates and non-local potentials, although feasible, is numerically 
demanding, and one therefore seeks simplified alternatives which nonetheless capture 
the essential physics; the remainder of this Section thus discusses such suitably- 
constructed alternative approaches. 

We start by discussing how these 'kinetic' contributions Ic and Inc lead to the 
introduction of collisional terms into the formalism: Firstly, one should obtain the 
exact equations of motion for both Ic and Inc by means of the effective hamiltonian 
of Eq. (114). Then, one should formally integrate these expressions, assuming their 
characteristic evolution in the basis specified by the unperturbed hamiltonian Hhfb 
can be described in terms of suitable dressed eigenenergies (here HFB quasiparticle 
ones); only at that stage should one impose any decoupling approximations. In other 
words, while the mean field approximations of Eqs. (51), (53) are used to define the 
unperturbed hamiltonian, in the final expressions one imposes generalized decoupling 
schemes (such as the one of Eq. (52)) to reduce correlations of the form {c^c^c^ccc) 
to pair operator averages, and thus obtain a closed system of equations in second 
order perturbation theory. Following the usual procedure implemented in the kinetic 
theory of dilute gases, one simultaneously simplifies the resulting expressions by 
assuming that collisions are well-separated in time, so that the mean fields evolve 
slowly compared to the collisional duration; the pseudopotential approximation is 
also made on the appropriately introduced effective interaction (i.e. consistent with 
restricting the anomalous averages to low- lying modes to avoid double-counting certain 
interaction effects), for details regarding this procedure and for the precise expressions 
of the generalized decoupling approximations the reader is referred to [111]. This 
approach yields a set of self-consistent equations for the variables (j), fh and h which 
include both damping mechanisms and collisional integrals. 

We are now ready to give our final second order theory in the effective interaction 
strength g. However, as Eqs. (116) have been explicitly formulated in terms of single- 
particle operators, the final expressions for the collisional terms given in [111] are quite 
lengthy. We have therefore chosen not to give these expressions here in full, but to 
focus instead on the simplest illustrative application of this formalism which highlights 
the essential physics of such contributions. 

Although one of the authors (NPP) was instrumental in the preceeding 
development (in collaboration with Burnett and Stoof) [77, 97, 88, 111], the final 
expressions for the full collisional integrals mentioned above were first given by 
Walser et al. in [110, 112]. This latter treatment is based on more 'conventional' 
formulations of a quantum kinetic theory, whereby one defines a set of suitable slowly- 
vaxying gaussian 'master variables' (corresponding precisely to the condensate and the 
normal and anomalous averages introduced earlier) whose evolution is subsequently 
studied. The full second order collisional integrals are given in Appendix A using the 
notation introduced in these latter works. Note that despite their differences in formal 
development and notation, the theory presented in the Appendix is identical to the 
presentation given in this section, as explicitly demonstrated in [111]; furthermore, 
both theories are formally equivalent [113] to the non-equilibrium Green's function 
approach originally put forward by Kadanoff and Baym [43] and subsequently applied 
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to trapped gases by Imamovic-Tomasovic and Griffin [114, 103, 104]. The full version of 
this theory has been apphcd (in the ergodic approximation) to the study of spherically 
symmetric harmonic traps [112], with appropriate generalizations introduced to study 
quasiparticle damping and finite coUisional duration [115], and recent extensions to 
the study of one-dimensional gases [116, 117]. 

In order to highlight the main physical issues of such theories, we now focus 
on a simple model. Following the discussion of Proukakis and Burnett, we recall 
that the exact equation of motion for the condensate mean field amplitude in level i, 
corresponding to Eq. (112) is given by [97] 

i^'^ = ^{n\ho\k)zk 

k 

* + 2pjiZk + {c\cjCh)^ . (117) 

ijk 

We also assume an idealized weakly-interacting Bose gas, for which (i) the conden- 
sate occupies only the lowest mode of the trap, denoted by the label '0', and (ii) 
the thermal particles are diagonal in the bare particle basis, i.e {n\ho\k) = SnSnk, 
where is the eigenenergy of trap level n. Although the theories of Walser et 
al. [110, 112] and Proukakis [111] explicitly contain the pair anomalous average 
m(r) — ipi{r)(pj{v)Kij in a self-consistent manner (see Appendix A) we simphfy 
our current presentation even further by assuming that (the static values of) such 
averages can be ignored. In this idealized limit we find the following contributions 
to the evolution of the condensate and the non-condensate (see, e.g. [Ill, 118, 119]), 
which are shown diagramatically in Fig. 4: 

Condensate Evolution: Beyond the usual HFB contributions (condensate mean 
field, normal and pair anomalous averages) of Eq. (112), we find the additional term 

^ Jill ^0ijk{c\cjCk) . (118) 

ijk 

This latter quantity, which has no contribution to the condensate mean field to first 
order in the interatomic potential, i.e. no static value within HFB, acquires a non- 
zero value in second order perturbation theory. In the context of the simplified model 
considered here, its respective equation of motion contains (after imposing suitable 
generalized decoupling approximations to reduce correlations of the form {c^c^cc) and 
{c'c^c'ccc) to simpler pair correlations (c^c)), among other terms, the contribution 
(see [77, 97, 88] for details) 

i^j^{c\cjCk) = {sj + Sk-ei) + --- 

+ 2V^feio [ni{nj + l)(nfe + 1) - (n, -|- l)njnk] zq. (119) 

We can now formally integrate this equation of motion. In doing so, we assume that the 
chosen mean fields do not vary appreciably on the timescale of a single collision; here we 
implicitly make an assumption, common in kinetic theories of gases, that collisions are 
short in duration and well- separated in time, and that we are only interested in 'long- 
time' evolution of the system, i.e. we can neglect the effect of any coherences present in 
the initial state of the system on the subsequent system dynamics. This enables us to 
approximate the quantity Zo{t') appearing in the integrand by Zo{t') = e'^o(*~* ^^'^zoit), 
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(and corresponding inverse processes) 



Figure 4. (colour online) Schematic of characteristic coUisional processes 
involving either the transfer of a thermal atom into the condensate (left), or 
collisions within the thermal cloud leading to redistribution among thermal energy 
levels (right). Characteristic rates for these processes are also given; note that 
the inverse processes are also allowed, with their combined effect corresponding 
to Eqs. (122)-(124). 



which is called the 'Markov' approximation. In lowest order perturbation theory, one 
thus obtains [88] 

X 2 [n,{nj + l)(nfc + 1) - {n, + l)njnk] it)zQ{t) . (120) 

Assuming no dependence on the initial state and introducing the shorthand notation 
£jkiQ — Ej+Efc— ej~£o we can evaluate the above integral (by introducing a convergence 
factor) to obtain 

r dt'e-^'^/^^'^"^"^''''^ = 7^Sie,k^o) + iv(—) , (121) 

J-oo \£]kiO/ 

where V denotes the principal part, and the 'loss' of memory effects has enabled us 
to extend the lower limit of integration to — cxd. The presence of both a real and an 
imaginary contribution in the expression for (clcjCk) implies that the appearance of 
this term in Eq. (118) leads both to a change in the amplitude of the condensate mean 
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field zq (arising from the ^-function term) and a change in the system's frequency 
(principal value term). The former contribution reads 

^ = • • • + ) ^0 XI l^oyfcl^ 5 {eoijk) 

X [{n, + l)njnk - n,{nj + l){nk + 1)] . (122) 

Here Voijfe describes a process leading to the creation of a condensate atom (level 0) by 
the collision of two thermal particles (levels j and k) , with the excess energy carried 
by a thermal particle in level i (see Fig. 4). Such a process leads to a 'production 
factor' which is proportional to (uj + l)njnk with its inverse process, also included in 
Eq. (122), yielding a corresponding factor of ni{nj + l)(nfc + 1). Finally, consistent 
elimination of the high- lying states (see e.g. [78, 119]) enables us to replace the exact 
interatomic potential V by an effective T-matrbc; in that case, the interaction strength 
appearing here can be treated as a contact potential of effective strength g (combined 
with a self-consistent high-energy truncation in the definition of the anomalous aver- 
age). 

Non-condensate Evolution: Following similar arguments, the evolution of the 
thermally excited population rii is found to contain (in addition to other terms) the 
following two contributions (shown schematically in Fig. 4): 

(i) The factor 

^ = '"+ [^jJ2^Voijkf5{eoijk) 

X |2;op [{n, + l)njnk - n,{nj + l)(nfe + 1)] (123) 

clearly corresponds to the same coUisional process as described earlier for the 
condensate, Eq. (122), and displays the effect of particle exchange collisions between 
the condensate and the thermal cloud on the population of a given single-particle 
thermal level. 

Observation of these equations shows that while the scattering of a particle into 
state i is bosonically enhanced by the factor (n^ + 1), the corresponding scattering 
into the condensate does not feature spontaneous growth; in other words, one obtains 
l^op both in Eq. (123) and in the evolution of the condensate population, d\zQ\^ jdt, 
arising from Eq. (122). This highlights the fact that the condensate mean field cannot 
grow from zero initial value, a well-known limitation of any mean field theory. (This 
issue can be cured within the context of stochastic approaches, briefly reviewed in 
Sees. 6.2.2-6.2.3.) 

(ii) The factor 

dn; f Ait\ 



dt 



^ ^ ^ \^rnijk\^ HSmijk) 



mjk 

X [{ni + l){nm + 'i-)njnk-ninm{nj + l){nk + l)] (124) 

arises from scattering of particles from states j and k into states i and rn, and the 
inverse process, and has associated with it the usual bosonic enhancement ((n. -|- 1)). 
This term is present irrespective of the existence of a condensate and corresponds 
to the usual scattering factors and amplitudes appearing in the coUisional integrals 
of the (classical) Boltzmann equation for scattering of particles in a thermal gas 
[120, 121, 122]. 
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Let us now revisit the above discussion in a more realistic context: Our preceding 
toy model has been expressed in terms of bare single-particle energies. In general, 
populations (c-Cj) are however not diagonal in this basis, resulting in a plethora 
of additional terms in the corresponding dynamical equations within such a basis 
(see Appendix A). We stress that the appearance of such contributions is formally 
correct, and any self-consistent treatment which includes them is not restricted to the 
extremely weakly-interacting regime gn <C fko; nonetheless, when one chooses such a 
basis, the interpretation of the resulting equations may not be as straightforward. 
In order to meaningfully interpret all contributions arising in this approach, one 
should ideally transform all expressions in terms of quasiparticle populations, which 
can be physically identified as the corresponding appropriate dressed particles due to 
interactions. Such a reformulation was performed (in the so-called Popov limit) by 
Imamovic-Tomasovic and Griffin in the context of non-equilibrium Green's functions 
[103], while the theory of Walser et al. [110] was also recast in more compact form 
using similar notation by Wachter in [123, 124]. 



4-4-^- Perturbative formulation beyond an appropriately generalized basis: In our 
mean field approximations (51), (53) which led to the HFB hamiltonian, we have only 
considered up to quadratic mean fields, i.e. averages of two non-condensate operators, 
namely fi = {S^S) and m = (SS); the unperturbed hamiltonian Hmf of Eq. (113) 
was thus limited to the HFB hamiltonian. However, the exact equations of motion 
for (p (Eq. (88)) and 6 (Eq. (90)) identify the crucial role of the triplet {S^S6) which 
was ignored in much of the early mean-field literature (but see also [48]). In order 
to also include this, one could thus consider resorting to a more general decoupling 
approximation for S^S^SS than the one given by Eq. (51). In particular, one could 
impose an approximation of the form [88, 111] 

6^^S ~ A{d^)6^d + {S^^)6S + {6S)S^^ - [2{6^){5U) + {66) {S^^)]. 

+ 2U6^6^5)5+ {6^6)6^] . (125) 

Following our previously introduced notation, substitution of this generalized 
mean-field approximation into the expression for H4 (Eq. (28)) would generate, in 
addition to the contributions {5Hi^ + 6H^'^'^) of Eq. (54), an extra contribution to 
the linear part of the hamiltonian J dT6'^{- ■ ■) + h.c, of the form 

^HTRip J ^j-s^s^sS) + h.c. . (126) 

We shall now follow our earlier perturbative prescription of Eq. (113) {H = 
Hmf + {H — Hmf))', however, this time we choose a slightly more general unperturbed 
hamiltonian which also includes 5Hi^^^ , i.e. we choose Hmf = Hhfbt where we 
have identified a new 'Hartree-Fock-Bogoliubov- Triplet' basis via (making explicit the 
split into Hartree-Fock and Bogoliubov contributions) 

Hhfbt = Hhfb + 6H™'' = {Ho + 5H^^ + SH^^^ + 5H™^) 



+ 



(127) 



The above expression additionally includes a mean field correction 5Hq^^^ = 
g J dr{6^66)(j)* + c.c. to the zeroth-order hamiltonian arising from an appropriately 
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generalized form of Eq. (53) which additionally includes (5^5S) in its right hand 
side. Since such mean field contributions merely shift the ground state energy of the 
system and have no additional effects on the dynamics of interest, we shall henceforth 
disregard in our discussion all non-operator mean field contributions SHq^, 5Hq'~"^ 



It is easy to see that the above expression still constitutes a generalized quadratic 
hamiltonian. Hence, this hamiltonian can in principle be diagonalized by means of the 
Bogoliubov transformation discussed earlier; the first such discussion was presented 
by one of us (NPP) in early work [88] in an attempt to justify from a microscopic basis 
the generalized HFB theory of Sec. 3.4, in which pair and triplet anomalous averages 
play a key role in upgrading the effective interaction between two atoms into a many- 
body one. However, it was found in [88] that such a theory cannot be microscopically 
derived by variational methods, as had been previously argued (on slightly different 
and more general grounds) by Bijlsma and Stoof [86]; the reason is that one runs 
into the same difficulties regarding the inconsistent treatment of atomic interactions, 
as in the full HFB theory. (This is in fact what led to the 'heuristic' addition of the 
generalized effective interaction .g(r).) One way to overcome this complication is based 
on a minor redefinition and handling of the hamiltonian; this is effectively what was 
done by Zaremba, Nikuni and Griffin, as discussed below - although their arguments 
were actually given in a slightly different manner. 

For ease of subsequent comparison, we give here the full expression for the above 
perturbing hamiltonian (ignoring energy shifts of the form SHq ), which takes the 
form 



4-4-3- Perturbative Distribution Function Formulation: Treatment of Zaremba- 
Nikuni- Griffin: The preceeding discussion demonstrated how use of perturbation 
theory facilitates a generalization beyond the static thermal cloud approximation, 
even after imposing the mean field approximations. One can actually use this 
approach to formulate equations which are similar in spirit to those encountered in the 
kinetic theory of classical gases, but are appropriately generalized below the transition 
temperature by the additional inclusion of the condensate mean field. 

Such an approach was formulated in the early 1980's in seminal work by 
Kirkpatrick and Dorfman [105, 106, 107], with related work undertaken by Eckern 
[108]. In particular, Kirkpatrick and Dorfman derived a closed kinetic equation for the 
quasiparticle distribution function of an inhomogeneous Bose gas below the transition 
temperature. The evolution of such a distribution function contains both 'streaming' 
and 'collisional' terms, and momentum enters explicitly as a system variable. This 
approach was recently employed and extended by Zaremba, Nikuni and Griffin [109] , 
as outlined below. In terms of the 'classifications' introduced earlier, this theory can 
be thought of as the consistent time-dependent extension of Hartree-Fock-Bogohubov- 
Popov of Sec. 3.3.1 which additionally includes collisions within the thermal cloud and 
particle-exchange collisions between condensate and thermal atoms. 

We start our discussion by reformulating the preceeding problem in terms of a 
suitable density matrix as follows: we assume that the system is initially described 



and 6H™P. 



H" = H — Hhfbt 




(128) 
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by a density matrix p{to), such that the expectation value of a general operator O is 
given by 

(6) = Tvp{U)d{t) = Tvp{t, to)d{to) . (129) 

Here Tr denotes the trace and we have introduced the density matrix p{t,to) = 
U{t,to)p{to)UHt,to) where C/ is a unitary operator defined by 



ih^lJ{t,to) 
at 



The corresponding density matrix evolves according to the equation of motion 
d 



H,ait),p{t,ta) 



(130) 



(131) 



and the problem thus reduces to the identification of the appropriate effective system 
hamiltonian /feff- 

In this section, we follow closely the discussion given by Zaremba, Nikuni and 
Griffin in [109]: In their work, they chose H^s such that the selected hamiltonian 
reproduces accurately the exact equations of motion for the condensate and the non- 
condensate operator, given respectively by Eqs. (88) and (90j. As before, the effective 
hamiltonian was split into an unperturbed hamiltonian (iJunp) and a contribution 



to be treated in second order perturbation theory, via Hes — -H^unp 
perturbing hamiltonian, H'", defined as [109] 

9 



H'", with the 



6'' {2h(j) + rh(j)*) + {2n<f)* + m*<l)) 5 



I j drS^^S -2g j 

-g j dv\^5'^{5'^55) + {5'^5'^5)5 



dv 



(132) 



It is instructive to compare this perturbing hamiltonian H'" to the one introduced 
earlier in Eq. (128). To enable a more direct comparison, we separate here the original 
H2 hamiltonian of Eq. (26) into two contributions (see also Table 1) , which we shall 
henceforth refer to as 'Hartree' and 'Bogoliubov' terms, via 



Ho = H. 



H 



jjBOG 



We find that H'" , is related to the perturbing hamiltonian H" of Sec. 4.4.2 via 



H" 



(5H, 



bog\ 



HF 
2 



5Hi 



TRIP\ 



H" 



H. 



BOG 



(134) 



The particular choice of the perturbing hamiltonian directly identifies the 
unperturbed hamiltonian which defines the spectrum of elementary excitations; here 
we clearly see that the main difference to the discussion of Sec. 4.4.2 is the contribution 

9 
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whose role is to change the excitation energies from single-particle to (finite 
temperature) quasiparticlc ones. 

If one were to literally follow the above argument, the unperturbed hamiltonian 
within the context of the Zaremba-Nikuni-Griffin approach should be given by 



-f^unp = H — H'" = Hq 



HF , XTJBOG , XTJTRIP\ 



Writing the single-operator contributions within ffunp explicitly as 



(136) 



/ 



drS^ (ho(t> + g\(t)\^(t> + 2gn(l) + gnuj)* + g{S^6)'^ + h.c. , (137) 



and observing the exact time-dependent equation for the condensate mean field given 

by Eq. (88), wc immediately see that these terms provide us with as expression for 
the time-dependent chemical potential. To make direct contact with the original 
presentation Zaremba, Nikuni and Griffin [109] , we use the Madelung transformation 
4>{r,t) = |0(r, f)|e*^('"'*'' to obtain the hydrodynamic equations (given in full in Sec. 
4.4.4) for the condensate density no = \(j){r,t)\'^ and Vc(r,i) = {h/m)V0{r,t). This 
procedure identifies a local chemical potential Hc, of the form [109] 

Aic = /io + 2,9n + j^Re |m(r )' + (^^M)^ } ; (138) 

here we see that that contributions arising from the anomalous averages provide a net 
real contribution to Hc, as physically relevant. These 'anomalous' contributions can 
change the chemical potential by introducnng modifications to the coUisional amplitude 
due to the medium in which the collisions are taking place (see discussion of mo(r) 
in Sec. 3.4). To avoid consistency issues with the unified treatment of interactions 
when simultaneously imposing the pseudopotcntial approximation (already implicit 
in this treatment), Zaremba, Nikuni and Griffin limited their analysis to the 'Popov' 
limit [67], in which all anomalous averages are assumed to be zero; this should not be 
misinterpreted as indicating that suc;h terms have no effect, as in fact an important part 
of their role is taken into account in introducing particle-exchange collisions into the 
theory. This procedure is essentially equivalent to maintaining interaction effects in the 
chemical potential and excitation energies only to first order in the effective interaction 
strength g, whereas the effect of interactions in calculating collision integrals is treated, 
as required, to second order in g [109]. In this limit, the condensate chemical potential 
acquires the form 

^2y2 

Mc = Mo + 'igh = — - — ^ + Fext + gino + 2n) . (139) 



We now also consider the quadratic hamiltonian {H2 + 5H2^) appearing in Eq. 
(136). It is easy to see that, as argued by Zaremba, Nikuni and Griffin [109], this is 
directly equivalent to the Hartree-Fock hamiltonian, namely 

In other words, the single-particle excitation energies are dressed by interactions in 
the manner discussed in Sec. 3.1.2. Alternatively, one can think of the condensate 
and thermal densities as providing an additional mean field potential through which 



2^ , [Fe,t(r,t) + 25(|<^|2 + n)] 



(140) 
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the atoms propagate. This leads naturally to the definition of a generalized mean field 
potential given by 

U{v,t) = V,^i{v,t)+2g[\cl>{v,t)\^ + h{v,t)\ . (141) 

Note that use of such a generalized mean field potential is common in theories based 
on Hartree-Fock eigenenergies; for the usual harmonic traps under consideration, and 
repulsive effective interactions, this takes the typical form of a double-well potential. 

Having identified our perturbing hamiltonian, we use Eqs. (129)-(131) to express 
by formal integration the expectation value of a general operator O as 

{d{t)) = Tvp{to) ^ul{tM)d{to)MtM-i f^dt'ulit'M) 

X [ul{t,t')6{ta)UaM , H"'{t')\ Uo{t',to)} ■ (142) 

In this expression, the first term depends on the initial correlations, while the second 
term can be identified as dynamical collisional effects arising from the perturbing 
hamiltonian H'". 

We are interested in using Eq. (142) to evaluate anomalous correlations of 
non-condensate operators, which we shall assume to be negligible initially, i.e. we 
can consistently ignore the first contribution. As the evolution of U is number 
conserving, in evaluating such anomalous averages one need only maintain terms 
which explicitly conserve the total atom number (thus discarding non-atom-number- 
conserving contributions arising here). We also assume as before that the collisional 
duration is short compared to characteristic evolution times of the parameters we 
are interested in. This enables us to make the Markov approximation Uo{t,to) ~ 
exp{-iHHF{t){t - to)}. 

As mentioned earlier, the novel feature of our current discussion is the formulation 
of the problem in terms of a distribution function /(p,r,t), for which the description 
of the gas can be seen as a direct extension of the kinetic theory treatment of a 
classical gas [125] below the transition point. In their original presentation [105, 106], 
Kirkpatrick and Dorfman formulated their discussion in terms of the full quasiparticle 
distribution function, although for most temperatures of experimental relevance for 
trapped Bose gases+ a discussion in terms of dressed single-particle eigenstates should 
suffice. Thus, the formulation of Zaremba, Nikuni and Griffin is based on a single- 
particle Wigner distribution function. 

For an atom of momentum p, at location r and time t the distribution function 
/(p, r, t) is defined as the expectation value 

/(p,r,t) = (/(p,r,i)) (143) 

where /(p, r, t) is the Wigner operator 

/(p, r,t) = J rfr'e^P-'/'^jt + ^ , to) <5 (^r - ^ , to) • (144) 

The above formulation enables us to obtain the non-condensate density via 

Hr,t) = J (2^/(P,r,t). (145) 

+ The usual criterion for this to be valid is that naX^g <C 1, which is certainly well-satisfied in 
current dilute atomic gas experiments. 
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Momentum variables are introduced into the system by expanding the non-condensate 
operators into Fourier modes (over non-zero momenta) [106, 107] 

5{v,t^) = J^^Cpe^P-, 5\vM) = -^Ect.e-P- . (146) 
p 



We now proceed explicitly to calculate the value of (5'^ 85). The temporal integrals 
arising in the second contribution of Eq. (142) are evaluated as in Eq. (121) via 

^t>eiiec+i2-i3-i4) p. ^^(^^ + £2 - £3 - £4) 



/ 

J — ( 



(£c + £2 — £3 — £4j 

Importantly, due to the choice of Eq. (140) as the unperturbed system hamiltonian, 
the energies = £i{v,t) = £{pi,r i,t) appearing in this and all subsequent expressions 
are evaluated 'semi-classically' in the Hartree-Fock limit via (see Eq. (62)) 

£i(r,f) = ^ + Kxt(r„t) + 2g[|<^(ri,t)|2 + n(ri,t)] . (148) 

Correspondingly, the condensate energy. Eg, is given by 

£c = + ^m?;^ , (149) 

where /i^ is the condensate chemical potential of Eq. (139) and v^. denotes the 
condensate speed, which should be determined self-consistently from the quantum- 
mechanical current 

_h_\ - (t>*V(t) 

2^/ Hi}' 



Mr,t)=^[^ V'\J (150) 



We proceed by replacing momentum sums (1/1^) by integrals / dp/(27rfi)'^ 
with the Kronecker delta symbols arising from the commutation relations replaced 
by Dirac delta functions, i.e. 5p o ^ (l/^)(27''fi)^'^(p)- To obtain a closed system of 
coupled equations, we decorrelate higher-order correlations in the fluctuation operators 
Cp by means of Wick's theorem, and express these decorrelated averages in terms of 
the Wigner distribution function via (c^j^Cp^) ~ 5p^ ,p^ f{pi,ri,t) . 

Thus, in direct analogy to Eqs. (120)-(121) [77, 88] the triplet contribution 
acquires to lowest order in g the value [109] 

X {2TThfS (mvc + P2 - Ps - P4) S {sc + £2 - £3 - £4) 

X [.f2(.f3 + l)(,f4 + 1) - (./2 + l),/3,/4] (151) 

where /j = /(pi, rj, t). Note that the same procedure can in principle be used to derive 
an expression for the pair anomalous average m(r), whose real contribution yields the 
correction to the scattering amplitude due to many-body effects discussed in Sec. 3.4. 
The equation for the condensate evolution (Eq. (88)) now becomes 

ih^{r,t)= [ho + g{\<j){r,t)\^ + 2n{r,t))-iR{r,t)'\^{r,t) (152) 
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where the contribution iR arising from the triplet correlation function describes the 
exchange of particles between the condensate and the thermal cloud (c.f. Eq. (122)) 

^, , {5'^U){Y,t) 

^2 2 /■ ^P2 [ dps I' rfp4 

J {2nh)^ J {2Trh)3 J (27r/i)3 

X {2TTh)^S {mVc + P2 - P3 - P4) <^ (£c + ^2 - h - ^4.) 

X [/2(/3 + l)(/4 + 1) - (/2 + l)/3/4] (153) 

The term iR can be interpreted as follows: Out of equilibrium this term leads 
to transfer of atoms between the two subsystems; once local equilibrium is (re- 
)established between the condensed and thermal components, iR = identically, and 
the condensate is described by the simpler Hartree-Fock expression of Eq. (96). 

Having identified the 'growth' term in the condensate evolution equation, we 
must also derive the corresponding evolution of the thermal cloud in order to obtain 
a closed system of equations. Following the usual formulation of the kinetic theory of 
gases [105, 106, 107, 108], the evolution of /(p, r, t) in the presence of a slowly-varying 
external field U{r, t) (for which we can make a gradient expansion) takes the form 

+V.Vr/-(VrC/)-(Vp/), (154) 

where Vr and Vp denote gradients in position and momentum variables respectively. 
Although in the absence of collisions / would obey Liouville's equation df/dt = 0, 
in the presence of collisions, the value of / changes at a rate C[f], which is called 
the collision integral. Thus, the problem of transport in a gas reduces to the 
identific;ation and calculation of the appropriate c;ollision integrals. The effect of 
collisions between the atoms is included into the evolution of the distribution function 
via the contribution 

f = \Trp{t, to) [/(p, r, to) , H"'{t)] . (155) 
dt in I i 

This yields the following expression for the evolution of the distribution function 

|[ + V • Vr/ - ( Vrt/) • (Vp/) = Ci2 [/] + C22 [f] ■ (156) 

Note that due to the presence of the mean fields of condensed and thermal 
atoms, all atoms experience the generalized potential U{r,t) defined by Eq. (141). 
Equation (156), first derived in slightly modified form in this context by Kirkpatrick 
and Dorfman [106], is known as the 'Quantum Boltzmann Equation' (see, e.g., 
[120, 121, 122]); this equation reduces to the classical 'Uehling-Uhlenbeck' equation 
[126] when Ci2[/] = 0, i.e. above the transition temperature when there is no 
condensate present. To strengthen the link to our introductory discussion (Sec. 4.4.1), 
we note that the coUisional integrals Ci2[/] and C22[/] appearing in Eq. (156) are the 
direct analogues of Eqs. (123) and (124) (and their more general expressions given in 
[111, 110] in terms of single-particle variables), now expressed explicitly in terms of /. 
They describe the following processes: 

• A collision involving the transfer of an atom from the thermal cloud into the 
condensate and its inverse process (c.f. Eq. (123)) 

2i j|2 f d.P'i f f dp4 



^ r^i ^'^ 2in2 f d.p2 f dp3 f 



CONTENTS 



52 



SELF-CONSISTENT TREATMENT OF COLLISIONAL INTERACTIONS 




NON 
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Non-condensate Evolution 



Quantum Boltzmann 
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+ , 



Figure 5. (colour online) Sciiematic of characteristic collisions included in the 
self-consistent Gross-Pitaevskii-Boltzmann formalism of Eqs. (152)-(156), often 
referred to as the 'ZNG' formalism . 



X {2'Khf5 {m^fc + P2 - P3 - P4) X (5 {Sc + 62 - £3 - £4) 

X {2TThf [5(p - P2) - 5(p - P3) - <5(P - P4)] 

X [(/2 + l)/3/4 - /2(/3 + l)(/4 + 1)] . (157) 

Clearly this term must be related to the condensate source term (c.f. Eq. (122)), 
and indeed one finds 

R{v,t) 



dp 



\cj^{r,tWJ (2^;.)3^^-— 

As a result, the final theory does actually preserve the total number of atoms in 
the system. 

• A collision between two thermal atoms (c.f. Eq. (124)) given by 



C22[/] 



47r 



-9 



dp2 



dps 



dp4 



£4) 



h" J {2T:hf J {2TThY J {2T:hY 
X {2-Knf5 (p + P2 - P3 - P4) (5 (£ + £2 - £3 

X [(/ + l)(/2 + l)/3/4 - //2(/3 + l)(/4 + 1)] . (159) 

We have thus arrived at a self-consistent second order description of a partially 
Bose-condensed gas in terms of a dissipative Gross-Pitaevskii Equation (Eq. (152)) 
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for the condensate mode, coupled to a dynamical reservoir of thermal atoms obeying 
a 'Quantum Boltzmann Equation' (Eq. (156)), as depicted diagrammatically in Fig. 
5. These equations should be solved simultaneously with Eqs. (145), (141), (153), 
(157), and (159) with the energies of the condensate and its excitations determined 
self-consistcntly from Eqs. (149) and (148). An important advantage of this theory 
is that it incorporates thermal cloud dynamics in a self-consistent manner. In the 
context of generalized mean fields and symmetry breaking, this is the most advanced 
model that is currently amenable to direct numerical simulations. 

The first key step towards such a derivation was made possible by the pioneering 
work of Kirkpatrick and Dorfman [105, 106, 107]; another crucial step in the context 
of a mean field formulation was played by the identification of the role of the triplet 
correlation (5'^ 55) in the description of condensate kinetics (highlighted by one of 
us (NPP) in collaboration with Burnett and Stoof [77, 97, 88]). These elements 
were beautifully tied in together in important work by Zarcmba, Nikuni and Griffin 
[109], whose arguments were presented (in slightly modified form) in this section. 
Over the past decade, these latter authors, along with their collaborators (including 
one of us - BJ) have used these coupled equations extensively to study the non- 
equilibrium behaviour of ultracold Bose gases (details given below). In these works, 
they introduced the shorthand notation 'ZNG' (after their initials), and for notational 
convenience we shall occasionally refer to the theory presented in this section as the 
'ZNG' theory. For completeness, we note that the theoretical number-conserving 
description formulated by Stoof [127, 128, 129] and Gardiner and ZoUer [130] presented 
in Sees. 6.2.2-6.2.3 include the above equations as a limiting case of their formalism. 

In terms of the physics contained in these equations, depending on the coUisional 
rate, or equivalently on the coUisional mean free path, one can identify two distinct 
regimes in the response of the system: In the coUisionless regime, which encompasses 
most of the current experiments with ultracold Bose gases, the effects of the mean fields 
dominate and collisions can be treated as perturbative corrections; here the full Wigner 
distribution function /(p,r,t) is needed to accurately describe the thermal cloud. In 
the opposite coUisional regime, strong interactions lead to rapid local equilibration, 
and the thermal cloud can be described by a few local macroscopic variables, as in 
ordinary fluid mechanics; in our current context, however, these variables for the 
thermal cloud are also coupled to corresponding ones for the condensate parameters, 
thus leading to a hydrodynamic description which can be visualized as an extension 
of Landau's two-fluid hydrodynamics [74] . The formalism presented above has been 
used extensively both in the coUisionless and in the coUisional regimes (also by the 
authors - mainly BJ), and some examples are briefly mentioned below. 

An important application of this theory performed by Bijlsma, Zaremba and 
Stoof [131] addressed the issue of condensate formation (see Sec. 7.2) by relying on 
the ergodic approximation, which assumes that equilibration is rapid for atoms with 
similar energies, thus enabling the phase-space variable fi = /(r^, p^, t) to be expressed 
only in terms of energy variables. These calculations were subsequently generalized 
to strongly non-equilibrium regimes (and the requirement for ergodicity lifted) by one 
of us (BJ) [132, 133, 134, 135, 136, 137], by representing the phase-space density by 
a collection of N discrete 'test particles', with collisions between them handled via an 
appropriate Monte Carlo sampling scheme, thereby extending earlier work [138] . 

Importantly, coupling a dissipative Gross-Pitaevskii equation to a Quantum 
Boltzmann equation as in the 'ZNG' approach enables an assessment of the relative 
importance of the various coUisional processes arising in ultracold gases, as first 
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2.30 




Figure 6. Role of coUisional proceses involving thermal atoms on the condensate 
dynamics following its initial excitation into a monopole, or breathing, mode. 
Left: (a) Frequency {ojm) and (b) damping rates (P) vs. reduced temperature 
in the absence of any collisions involving thermal atoms {C12 = C22 = 0; 
filled circles), including thermal-thermal collisions (C12 = 0, C22 7^ 0; open 
circles) and including all coUisional processes (C12, C22 7^ 0; triangles), where 
u>o is the radial trap frequency and is the critical temperature of the non- 
interacting gas. Right: Corresponding damped oscillations of the condensate 
width at T/T^ = 0.64 in the absence of any thermal collisions (dashed) and with 
all coUisional processes included (solid), corresponding respectively to the filled 
circles and triangles in the image on the left. Results are based on 5 X lO" **^Rb 
atoms in an isotropic trap of frequency tJo = 2tt X 187 HZ (Reprinted figures with 
permission from B. Jackson and E. Zaremba, Pliys. Rev. A 66, 033606 (2002). 
Copyright (2002) by the American Physical Society.) 



discussed by Eckern [108], and analyzed more extensively by Zaremba, Nikuni and 
Griffin and their collaborators (including one of us - BJ). These coUisional processes 

involve both particle exchange collisions between the condensate and the thermal 
cloud, and 'redistribution' collisions within the thermal cloud. A typical example 
of their effect on the excitation of a 'breathing' mode, in which the condensate 
exhibits damped radial oscillations, is shown in Fig. 6; the observed increase in the 
damping rate with temperature should be contrasted to the undamped oscillations 
predicted by the GPE. Although the above formulation ignores many-body effects 
encompassed in the static value of the anomalous average (since rho = here), it 
nonetheless yields remarkable agreement with various experiments, as evident from 
the study of the dynamics of various types of excitations, such as scissor's mode [132], 
quadrupole excitations (see Sec. 7.1) and transverse breathing modes of elongated 
condensates [134]. The 'ZNG' theory also made the first quantitative predictions of 
vortex nucleation at finite temperatures [139]. 

Moreover, this theory was recently used by the authors to study the finite 
temperature dynamics of dark solitons [140] and vortices [141]. In the Hannover 
experiment [142], performed at a temperature of « 0.5Tc the phase imprinted soliton 
was found to decay when it reached the edge of the thermal cloud, a result perfectly 
supported by our numerical simulations (see Fig. 7). Application of this formalism to 
vortices shows that the core of an off-centred vortex spiralling out of the condensate 
due to dissipation is filled up by thermal atoms which closely follow the vortex 
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Figure 7. Decay of a dark soliton due to finite temperatures, as witnessed in 
a sequence of expansion images, corresponding to a different soliton propagation 
time in the original trap (as indicated) before that atomic gas is released from the 
trap. The self-consistent Gross-Pitaevskii-Boltzmann, or 'ZNG' theory correctly 
reproduces the decay observed in the Hannover experiment [142] as soon as the 
soliton reaches the condensate edge (top images), whereas the same soliton in 
a T = condensate would execute undamped oscillations in the trap (bottom 
images). (Reprinted figure with permission from B. Jackson, N.P. Proukakis and 
C.F. Barenghi, Phys. Rev. A 75 051601(R) (2007). Copyright (2007) by the 
American Physical Society.) 



trajectory, highlighting the importance of a fully-dynamical self-consistent calculation 
[141]. 



4-4-4- Hydrodynamic Description at Finite Temperatures: We now give the 
corresponding finite temperature hydrodynamic equations. In particular, the 

continuity equation (cf. Eq. (33)) within the 'ZNG' formalism acquires a 'source 
term' from the imaginary part of the triplet contribution, i.c 



dno , , 2g 
^+V-(noVe) = -Im 



whereas the corresponding imaginary contribution of the pair anomalous average 
vanishes. This is coupled to a 'generalized Euler equation', given here by (cf. Eq. 
(34)) 

d_ 

di 



m 



(161) 



with the condensate chemical potential /i^ already defined in Eq. (139); note that, 
within the Popov approximation which is imposed in the 'ZNG' approach, the latter 
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chemical potential ignores the real contribution from the pair anomalous average which 
is responsible for certain many-body corrections. 

These equations have been used extensively to discuss hydrodynamic features 
of interacting Bose gases, and have most notably led to the derivation of Landau's 
two-fluid equations [143], and the Landau-Khalatnikov two-fluid equations with 
hydrodynamic damping, including explicit expressions for all the transport coeSicients 
[144]. 

4.5. Brief Summary 

Particle- exchange collisions between the condensate and the thermal cloud, as well 
as binary collisions within the thermal cloud, can be included in a generalized time- 
dependent mean field by going beyond the mean- field approximations of Eqs. (51), (53). 
This can be achieved by treating the remaining part of the full system hamiltonian, 
namely {H3 — 6H1) and {H4 — {SHq -\- SH2)) perturbatively, where 5Hi are the 
contributions arising from the generalized mean-field approximations. Various such 
formulations have been discussed, differing essentially in the partitioning of the system 
hamiltonian into an unperturbed part which specifies the nature of the excitation 
spectrum, and a perturbing hamiltonian which governs the system dynamics. 

Kirkpatrick and Dorfman were the first to obtain a Quantum Boltzmann Equation 
for the quasiparticle distribution function in the presence of a condensate. Combining 
this with the crucial inclusion of the anomalous triplet correlation (S'SS) (whose 
importance was also identified earlier by one of us - NPP - in collaboration with 
Burnett and Stoof), Zaremba, Nikuni and Griffin successfully formulated an important 
^self-consistent Gross-Pitaevskii-Boltzmann' approach which they termed the 'ZNC 
theory. In this theory, which has been formulated such that the effective hamiltonian 
fully reproduces the exact equations of motion for the condensate wavefunction and 
the non-condensate operator, excitations are specified by Hartree-Fock energies. The 
resulting theory (which also arises as a limiting case of formalism yet to be presented), 
fully accounts for the dynamics both of the condensate (via a dissipative Gross- 
Pitaevskii equation) and of the thermal cloud (via a Quantum Boltzmann equation) 
in a self-consistent manner; the predictions of this theory appear to be in excellent 
agreement with experiments both in the collisional and collisionless regimes. Despite 
ignoring many-body effects contained within the static anomalous average contribution 
rh§{r) (implicitly included in alternative formulations), this 'ZNG' theory constitutes 
the best numerically implemented mean field theory to date. 

5. Pheise Fluctuations cind Number-Conserving Approaches 

In our treatment thus far we have replaced the c;ondensate operator (}()(r, t) by the 
classical field <t>{v,t). As a result, the condensate is no longer invariant under global 
phase changes, and thus breaks this symmetry of the full system hamiltonian by 
acquiring a definite phase. Two important implications of this arc: (i) a proper 
inclusion of fluctuations in the condensate phase is precluded from the outset; and 
(ii) the total particle number in the system is no longer conserved. In this section we 
systematically address both these issues. 
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5.1. Low Dimensional Systems: Fluctuations in the Condensate Phase 

Under conditions of reduced effective dimensionality, Bose-Einstein condensation 
cannot occur in a honiogcncous system [50, 51]; instead, quantum degeneracy leads 
to the appearance of a so-called 'quasi-condensate', which can be thought of as a 
condensate with a fluctuating phase [49]. Such low-dimensional geometries can now 
be routinely produced in numerous laboratories, under c;onditions of tight confinement 
in one or two directions, which effectively freezes the motion of the atoms in the 
remaining tightly-confining direction(s) (see, e.g. [145, 146, 147]). We stress that our 
present discussion is limited to weakly-interacting systems whose kinematic behaviour 
is low-dimensional, but scattering can still be described as effectively three-dimensional 
[148], i.e. we are still far from the strongly-correlated 'Tonks-Girardeau' regime arising 
in purely ID Bose gases [149] (whose treatment lies beyond the scope of this Tutorial) 

For such systems, one can define two characteristic temperatures [149]: a 
'degeneracy' temperature signalling the onset of macroscopic occupation of a given 
quantum state (i.e. the low- dimensional analogue of the critical temperature); and 
a 'phase fluctuation' temperature signalling the transition to a regime where the 
phase of the 'condensate' begins to fluctuate. Our preceeding discussion is based on 
symmetry-breaking, which automatically assigns a phase to the condensate; as a result, 
fluctuations in the phase cannot be treated properly, and so far our analysis has only 
accounted for the local gradient of the condensate phase (superfluid velocity) ; in order 
to account for phase fluctuations properly, as required for systems of dimensionality 
less than three, a modification of the previous formalism is required. Note that phase 
fluctuations may also arise in three-dimensional geometries in the early stages of the 
formation of a Bose-Einstein condensate, although their effect is only expected to be 
important in a very narrow temperature region near the critical point. 

Fortunately, the inclusion of phase fluctuations can be done fairly straightfor- 
wardly, essentially by replacing Eq. (20) by the more general expression [150] 

*(r, t) = v/^^;^(M)e^^"('-'*) + V>'(r, t) . (162) 

To allow for the fact that the 'condensed' part of the system (first contribution) does 
not have a unique phase, we have introduced here an operator ^(r, t) to account for 
fluctuations in the phase of the 'coherent' part of the system. This quantity is useful 
when looking at correlations of such operators, e.g. evaluating correlations of the phase 
(contained in the one-body density matrix) at different locations, or times; thus, Eq. 
(162) provides a mathematical tool for computing correlation functions of the gas. 
Density fluctuations have also been included by introducing a corresponding operator 
tl>'{r,t). The low temperature Hmit of this equation was discussed in [149], with the 
current extension presented in [150], with one of us (NPP) playing an instrumental 
part in its numerical implementation [151, 152, 153]. As this issue digresses slightly 
from the main emphasis of this Tutorial, we only give some general arguments here, 
and refer the reader to the above-mentioned works for more details. 

To explain how Eq. (1G2) leads to the inclusion of phase fluctuations, we disc;uss 
how it relates to the symmetry-breaking expression of Eq. (22): In our preceeding 
analysis, we have chosen to describe the system in terms of correlations of fluctuations 
about the condensate mean field up to quadratic order, i.e. ignoring (or only treating 
perturbatively) the contributions from and H4. This is equivalent to making a 
lowest order expansion of the exponential in Eq. (162), such that e'^^""'*^ Ri 1-f i9{r, t), 



CONTENTS 



58 



In this limit, Eq. (162) can be re-cast in the more famiUar form 
*(r, t) = \/no(r,i) + [i x/no(r, i)^(r, t) + ^/.'(r, t) 

= (/)(r,i) + (5(r,t) . (163) 

Why is such an approximation insufficient for describing low-dimensional systems? 
Based on Eq. (162), the total system density at position r and time t becomes 

(*t(r,t)^-(r,i)) = no(e-'^>'*)e'^>'*)) + ([V>'(r,t)]tvi'(r,i)) 

+ ^/no(e-'^(■■■*)V''(r,^)) + ^([^'(r, i)]te'^"(--'*)). (164) 



The top line should reproduce the total system density, with the bottom line 
vanishing under the assumption of no correlations between the condensate and 
the non-condensate. Indeed locally this yields the identity n{r,t) — no{r,t) + 
h{r,t). If however one were to make the approximation e'^^'''*^ 1 + i9{r,t), 
then the system density would incorrectly acquire an additional contribution from 
no{r,t){9{r,t)9{r,t)). Although the latter contribution is small for three-dimensional 
homogeneous condensates, it blows up as p — > for dimensions d <2 (except at T = 
and d = 2), thus leading to well-known infrared divergences in the expression for the 
total density [49, 50, 51]. 

The above discussion (which has been put on firmer ground [151]) suggests 
that one can proceed with the usual established (one-loop) mean field treatment 
[86], and nonetheless obtain the correct equation of state by omitting the infrared- 
divergent contributions; this is based on the physical argument that such contributions 
have only entered the expression for the total density as a result of inappropriate 
handling of the condensate phase fluctuations (which leads to double-counting) . This 
procedure has been shown to reproduce all known results for (weakly-interacting) 
Bose gases in one, two and three dimensions [151], thus providing a modified mean- 
field theory and a general equation of state valid in all dimensions. In fact, this simple 
trick enables the direct calculation of equilibrium densities in the degenerate regime 
[151], while the identification of the expression for {9{r,t)9{r,t)) enables the direct 
computation of (non-local) phase fluctuations, and thus of correlation functions. This 
is illustrated in Fig. 8 which plots both density profiles and correlation functions for a 
one-dimensional Bose gas at temperatures where phase fluctuations are present. Note 
that the predictions of this theory have been further shown to compare favourably to 
more exact stochastic treatments discussed in Sec. 6.2.2. 



5.2. Number-Conserving Approaches 

In this section we critically revisit the notion of Bose broken symmetry, and essentially 
repeat the analysis of Sees. 3-4 without ever resorting to such an approximation. 
This is important from a fundamental point of view, and leads to certain non-local 
corrections whose importance is still under investigation. 

The concept of symmetry-breaking is familiar in many branches of physics, 
including condensed matter, particle physics and optics. However, in the case of 
'closed systems', as is relevant for finite trapped atomic gases, and because atoms, 
unlike photons, can neither be created nor destroyed, the notion of broken symmetry 
becomes somewhat 'ill-defined'. This docs not mean that approaches based on 
symmetry-breaking are incorrect, and indeed such approaches turn out to be very 
useful, providing in many cases excellent agreement with experiments. However, one 



CONTENTS 



59 




Figure 8. (colour online) Left: Ab initio determined density profiles for a trapped 
one-dimensional Bose gas via the modified method proposed in this section (split 
blue/red profiles) and the stochastic method of Sec. 6.2.2 (continuous brown 
curves) for two different temperatures. Inset shows how a separation into 'quasi- 
condensate' (upper green) and 'thermal' (lower green) profiles within the inner 
'degenerate' region (z < Al^ for T = 500nk) enables the determination of the total 
density in that region (blue). Outside the 'degenerate' regime (red curves) the 
density is computed by a semi-classical equation of state for a thermal gas in the 
many-body T-matrix approximation; the small discontinuity in the plotted profiles 
arises due to the use of two different equations of state. (Densities plotted as gn{z) 
scaled to the typical harmonic oscillator energy huJz, where g is the effective 
one-dimensional coupling constant, while positions are scaled to the harmonic 
oscillator length 1^). Right: Normalized first order correlation functions g(^'(0, z) 
plotted against z with distances scaled to the corresponding quasi-condensate size 
for each temperature (often referred to as the 'temperature-dependent Thomas- 
Fermi radius' [151]). Displayed temperatures, from top to bottom: T = 10 (black), 
50 (red), 150 (green) and 500 nK (blue), all below the transition temperature. 
The loss of phase coherence characteristic of low-dimensional systems is evident 
in the decay of the correlation function to zero, even at distances smaller than 
the effective quasi-condensate spatial extent. (Simulations performed for 20000 
^•^Na atoms with uj^ = 2tt X 3.5 Hz - see [150, 151, 153] for further details). 



may not be entirely satisfied with this picture, based on the argument that at any 
one time the atomic state consists of a definite number of atoms which one should 
be able, at least in principle, to determine experimentally. For a closed system 
with a fixed number, N, of bosonic atoms, the average (A^|^|7V) — {N\^^N) — 
identically (since the bosonic field operators, change the total number of atoms 

in the system from N to {N ± 1), and {N\N ±1) = in the orthogonal number 
state basis). It is then conceptually hard to interpret a non- vanishing average of 
the field operators. Usually one overcomes this problem by representing the state of 
the quantum system by a coherent superposition of states with different numbers of 
particles, •••|A^— 1), \N), | -1-1) i.e. describing the system by a so-called 'coherent 
state' |a); such a model has proven very successful, for example, in describing the 
laser [154, 155, 156]. Although experimental knowledge of the atom number is usually 
obtained statistically (in the sense of ensemble averages from multiple realizations) it 
is still hard to see how shot-to-shot number coherences would build up [161]. The issue 
of whether broken symmetry is a 'necessity', a 'reality', or a 'convenient mathematical 
tool' (or none of the above) remains however somewhat controversial; for example. 
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recent work by Yukalov [70] concluded that spontaneous gauge symmetry breaking is 
the necessary and sufficient condition for the existence of BEC in any kind of system, 
based on the notion of representative statistical ensembles [159]. 

In any case, significant progress can in fact be made without invoking the 
concept of Bose broken symmetry, i.e. without relaxing the requirement of number 
conservation. It turns out that the physically more intuitive (yet mathematically more 
cumbersome) number-conserving approaches lead to essentially the same equations, 
but with some (non-local) corrections which are not always important. 



5.2.1. Essential Methodology: The problem considered thus far can actually 
be reformulated in terms of operators which have very similar properties to 
operators considered in symmetry-breaking treatments, but yield difl:erent physical 
interpretations. For these operators, the total number of atoms, N, in the system 
is explicitly conserved, whereas the number of excitations, U]^ is not. In such 
treatments, the number of condensate atoms is actually not a separate variable, but 
is determined self-consistently from Nq = N — J2k ^k- The first discussion of such 
number-conserving approaches was given by Girardeau and Arnowitt in 1958 [160], 
with this approach 'rediscovered' in a slightly modified form by C.W. Gardiner [157] 
and subsequently extended by Castin and Dum [158] and Morgan and S.A. Gardiner 
[161]. While all these approaches are based on the same underlying theme, their precise 
implementation of number-conservation is not identical. In our subsequent discussion 
we briefly review the main ideas behind these works, with particular emphasis given 
to the recent dynamical approach of Morgan and S.A. Gardiner [161]. 



5.2.2. Approach of C.W. Cardiner: This approach was based on defining an 
annihilation operator A for the total number of particles and a separate operator 
for phonons. In particular, the operator a^. which annihilates an excitation in mode 
k was defined in terms of the usual single-particle operators dfe by 

Oik = —F^oiak ■ (165) 



While A and thus N act on the total number of particles, the operators a.k 
actually commute with A^. This enables approximations to be made in the 
treatment of the phonon operators, without the need for a change in the total 
particle number. Although such a treatment leads to explicit number conservation, 
the newly-introduced phonon operators only satisfy bosonic commutation relations 
approximately (for large N) via 

4fe' - ^afeo-i' ~ 4fe' • (166) 



Oik,a\, 



N 

By an appropriate reformulation of the system hamiltonian, the Bose field 
operator ^(r, f) is expanded in terms of these new operators, with the hamiltonian 
split into different contributions according to the number of phonon operators present - 
analogously to the separation in terms of the number of fiuctuation operators present in 
Eqs. (24)-(28). Next, Gardiner introduced a redefined interaction parameter = gN, 
which provides a measure of the interaction energy of the system. Keeping the 
parameter gj^ fixed, while increasing the atom number N essentially corresponds to 
taking the thermodynamic limit of the system under consideration [158, 161]. As a 
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result, the separation of the full system hamiltonian into these contributions can now 
be visualized as a systematic expansion in powers of (l/y/N), via the expansion 

m=0 ' 

where calligraphic notation has been used to highlight the fact that the hamiltonians 

appearing in the above expression are given in terms of the redefined phonon operators 
of Eq. (165). Each of these hamiltonians Hm is now expressed in terms of g^, with 
the index m indicating the number of phonon operators present in the system, e.g. 

Ho = J dr[4>''{r)]*h^''{r) + J dr |<?i^(r)|* , (168) 

where the condensate wavefunction has acquired a superscript N to highlight that all 
calculations are performed in a number-conserving manner. 

As iV ^ oo, higher than quadratic terms become negligible, thus yielding a 
number-conserving hamiltonian which can be routinely diagonalized by a Bogoliubov 
transformation. Keeping only the first contribution leads to the time-independent 
GPE 

hocl>^{r)+g^ |</.^(r)|'<^^(r) = A0o^(r) (169) 

The parameter A appearing here is essentially the chemical potential of the system, 
only in this treatment it arises naturally as a Lagrange multiplier, and should therefore 
be computed self-consistently. 

This generalizes straightforwardly to the time-dependent case, with a slight 
variant of this approach presented in more detail in Sec. 5.3. 



5.2.3. Number-Conserving Static Finite Temperature Bogoliubov Equations: A 
related analysis, based on a slightly modified version of the number-conserving 
formalism of Sec. 5.2.2, was performed by Morgan [78] in the context of a self-consistent 
second order perturbation theory, which focuses on the change in the system energy 
when a single quasiparticle is added to a particular mode. The aim of this analysis 
was to determine the modified form of the static Bogoliubov equations under explicit 
number-conservation upon additionally including the back-action of the thermal cloud 
on the static properties of the condensate. 

In brief, this approach is based on the following reasoning: The entire system 
hamiltonian, Eq. (17) is re-expressed in terms of explicitly number-conserving 
operators, and broken down to contributions with different numbers of non-condensate 
operators in each term (c.f. Eq. (167)). In this treatment, the number-conserving 
operators are defined somewhat differently to Sec. 5.2.2, by au = 0ldk, where 
Po = {Nq + l)^/^ao and Nq = alcto is the operator for the condensate number. 
Similarly to Eq. (167) this leads to an expansion of the full system hamiltonian in 
terms of c;ontributions Tim expressed in order of decreasing factors of condensate 
atom numbers. In the limit of large Nq, such an expansion justifies treating the 
quadratic part of the hamiltonian exactly, while maintaining additional cubic and 
quartic contributions perturbatively, analogously to the symmetry-breaking beyond- 
HFB treatment presented in Sec. 4.4. In particular, one calculates the change induced 
in the system energy due to the effect of (Tils -|- 7^4) in second order perturbation theory. 
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Apart from the inclusion of normal and anomalous averages (and therefore many- 
body effects), this approach additionally includes effects arising from quasiparticle 
collisions in a self-consistent manner. The perturbative inclusion of Tis -|- 71^ leads to 
a change in the static properties of the system (e.g. condensate shape, energy), which 
induce corrections to the basic hamiltonian (Tio + 'Hi + 'H2), and such corrections are 
also treated perturbatively in this approach. Combined with the enforcement of total 
number conservation, this leads to the inclusion of finite size effects into the formalism. 
We wish to focus here on the main modifications induced by the requirement of explicit 
number conservation, and therefore omit the technical details expounded in [162]. 

The main conclusion of this work was the derivation of a generalized set of 
Bogoliubov equations of the form 

( >C^(r,ei) 
1 _[^iV(r,-e,)l* 



Al^(r,ei) 
>(r,-ei) 




(170) 



In comparison to the expressions encountered earlier, the operators and M.^ 
introduced here (denoted by L and M in the literature) contain additional terms. In 
particular we find 

£^(r, e,) = K-\ + 2g|0^|2 + 2gno{v) + Sl{r, a) 

«£(r) + (5£(r,e,) (171) 



M'^ir, a) = g[(b^{r)r + gnii{r) + dM{r, e,) 



(172) 



where the approximate sign has been introduced to highlight that the parameter A 
appearing in number-conserving approaches is a self-consistently determined Lagrange 
multiplier, rather than the chemical potential, fi. This is obtained from 

ho + 9 \ct>^{r)\' + 2gho{r)] {r) + gm^ (,^^(r))* = A<^o^(r) . (173) 

The terms 5C and 5M. are corrections due to the change in energy when a quasiparticle 
is created in a mode i. In addition to describing the interaction of two quasiparticles 
(Landau and Beliaev damping discussed in Sec. 4.3.1), such corrections also include 
processes referring to the simultaneous annihilation or creation of three quasiparticles. 
The latter are expressed in terms of contributions of the form [78] 



X! [Am(ei)(l + nk + rim) + Bkm{ei){nm - nk)] 



(174) 



where Akm{u) and Bkr, 



(e^) depend on the energies of such excitations. 



5.3. Dynamical Finite Temperature Bogoliubov Equations 

Wc now wish to construct a dynamical number-conserving approach from first 
principles, essentially by generalizing the arguments presented in the preceeding 
sections. We thus seek to construct a formalism in which the Bose- field operator ^{r,t) 
is expanded in a symmetry- preserving manner. Assuming this could be achieved 
in terms of a slightly modified non-condensate operator 6'^{r,t), we would like this 
operator ideally to have the following properties: (i) explicitly ensure orthogonality. 
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(ii) have zero average, i.e. {5^ ) = 0, (iii) exactly satisfy bosonic commutation relations 
[(j-'^(r, f)), 5^(r', = (5(r — r'), and (iv) guarantee conservation of the system total 
atom number. Unfortunately, not all requirements are mutually compatible, and in 
defining (5^(r,t) one is typically faced with a choice of relaxing one (or more) of 
the above conditions. The choice is often motivated by the intended use of the 
resulting formalism, and leads to a range of similar number-conserving approaches 
[78, 160, 157, 161, 158]. Once the most suitable form of 5^(r,i) has been identified, 
one can perform an expansion of the system hamiltonian in ascending powers of this 
new fluctuation operator, in direct analogy to Eqs. (24)-(28) and (167). 

We now discuss how such a number-conserving formalism can be constructed. 
Firstly, we wish to make a comment about orthogonality: In our symmetry-breaking 
discussion of Sees. 2-4, we have argued that the quasiparticle amplitudes ■Ui(r) and 
i)i(r) obtained from the solutions of the symmetry-breaking Bogoliubov equations 
(Eq. (40) and finite temperature generalizations) are orthogonal to the condensate; 
this statement was actually based on the relation / dv\^Q{v)ui(v) -\- (/)o(r)t;j(r)] = 0. 
However, Morgan has argued [78] that this is only true in a generalized sense, whereas 
each of these two integrals is not separately zero. This implies that the quasiparticle 
functions defined by the symmetry-breaking Bogoliubov equations are not individually 
orthogonal to the condensate (j)o{^)- (This is actually a direct consequence of the fact 
that the Bogoliubov equations are not strictly speaking necessary conditions for the 
hamiltonian to be diagonalized [163]). 

We now wish to develop a formalism which explicitly guarantees the orthogonality 
by construction. In doing so, we shall maintain the condensate operator oq in the 
condensate part of the Bose field operator defined by Eqs. (20)-(21). Orthogonality 
between condensate and non-condensate can be guaranteed by explicitly introducing 
projectors onto the two orthogonal subspaces. In particular, we define the condensate 
annihilation operator, ao, by the projection P of the full Bose field operator ^{r,t) 
onto the 'condensate state' (f)^{r,t), i.e. 

ao(t) = P*(r,t) = J dr(b''ir,t)^{r,t) , (175) 

and correspondingly we define a non-condensate operator by the orthogonal projection 

S{r, t) = Q*(r, t) = J dr'Q{r, r', i)*(r', t) (176) 

where g(r, r', t) = 5(r - r') - (/)^(r, t)[4>^ (r' ,t)]* . 

A key role in our subsequent development is played by the Penrose-Onsager 
criterion for Bose-Einstein Condensation [34]. This states that for a system 
described by a single-body density matrix p{r, r', t) = (^^(r', t)^{r, t)) the condensate 
wavefunction corresponds to the mode of the system which has the largest eigenvalue 
(which is macroscopically large). The condensate number, Nq, is then related to the 
above density matrix via 

7Vo0^(r, t)^J dr'p{r, r', t)0^(r', t) . (177) 

Next, we expand the right-hand side of this expression into an appropriately defined 
orthogonal set {^q (r),(^j^(r)} via the following decomposition 

*(r, t) = ao(t>^{r, t) + S{r, t) = aocp^{r, t) + J2 W • (178) 
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This generates four contributions, as given by Eqs. (179a)-(179(i); each of these terms 

can be reduced to a simpler expression as indicated by the arrows, and the arguments 
for such reduction are presented below [164]. We thus have 

j dr'[0^(r',t)]*0^(r',t)(4ao)</.^(r,t) ^ (aJao)<^^(r, i) (179a) 
+ Jdr'[<p''{r',t)rcl>''{r',t){al5{r,t)) ^ {al6{r,t)) (1796) 
dr'{S\r',t)ao)(l)^{r,t)(j)^{r',t) ^0 (179c) 



+ J dr'{S\r',t)6{r,t))(l)^{r',t) ^0. (179rf) 

Let us now explain how we have obtained the indicated simplifications: Firstly, 
the expressions appearing on the right hand side of Eqs. (179a)-(179&) are 
obtained directly by using the completeness relation J dr' [(j)'^ {r' , t)]* (p^ {r' , t) = 1. 
Moreover, expanding the non-condensate operator 5^{r',t) as X^j^o "^i ('''' 
contribution of Eq. (179 c) becomes J dr' J2^^oV^i^' ^t)^!'^ i'^' {T,t){dlao). This 
quantity is identically zero, because by construction (fi and cj)^ are orthogonal, i.e. 
/ dr'(p*{r',t)(j)'^{r',t) = 0. The same argument holds for Eq. {179 d), which is re- 
expressed as / rfr'^^^Q<^*(r',t)(^^(r',f)(at^(r,t)). With these simplifications, we 
arrive at 

7Vo0^(r, t) = (a|,ao)0^(r, t) + {dlS{r, t)) . (180) 

If we now multiply this by [4>^{r,t)]* and integrate over r, we find that Nq = {dido) 
because the other contribution identically vanishes due to orthogonality, since 

dr[,/.^(r,i)]*(aS^(r,t)) = ^ / dr[,^^(r, f)]Vi(r, t)(aSai) = .(181) 

The identification A^o = {aldo), when substituted back into Eq. (180) yields directly 
the result 

{dl5{v,t))=Q, (182) 

i.e. there are no direct coherences between the condensate and the non-condensate 
within our specified number-conserving scheme. Eq. (182) is a rather significant 
result: it tells us that we could potentially choose the quantity al^5{v,t) as our 
fluctuation operator; this is a plausible choice due to the fact that it has a zero average 
(aj(5(r,t)) = 0, directly as a consequence of the implemented orthogonality between 
the condensate and non-condensate subspaces and the Penrose-Onsager criterion for 
BEC [158, 161]. 

However, for our separation into condensate and non-condensate contributions 
to be most useful, there is one additional property that we would like the desired 
non-condensate operator 5^{T,t) to satisfy: such an operator should ideally scale (at 
least approximately) as the square root of the number of non-condensate atoms - this 
is effectively the same scaling as in the symmetry-breaking treatments. 

Recalling that, to lowest order dj ~ \/Nq (Bogoliubov substitution), we thus 
define 
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as the fluctuation operator in the number-conserving formalism presented here, which 
clearly satisfies {S^{r,t)) = 0. So, what have wc actually achieved by the above 
considerations? For example, in the symmetry-breaking picture, we also had a 
fluctuation operator with a vanishing average (6) = 0. The important additional 
feature in our current treatment is that the condition ((5^) = docs not arise in an 
ad hoc manner (as was the case up to now), but rather it manifests itself as a direct 
consequence of the orthogonality between the two parts of the system. 

We should also make a comment about our particular choice for the number- 
conserving fluctuation operator. Deflning the operator via Eq. (183) leads to only 
approximately-satisfled bosonic commutation relations via (c.f. Sec. 5.2.2) 

5^(r,t),(5^(r',i))t] = ^Q(r,r',i) - ^{S'' (r' ,t))^ S{r,t) 

«Q(r,r',i). (184) 
Different choices of ^^(r, t) have been discussed in the literature, e.g. (aj/ \/iVo)(5(r, t) 



-^0/ 

[158], or {al/\/^)d{r,t) [157] which is the position representation of the phonon fleld 
operator ctk of Eq. (165); each of these approaches has distinct beneflts and drawbacks, 
and a more detailed discussion of the main factors affecting the optimal choice are given 
in [161]. 

Having identified the fiuctuation operator 6^(r,t) in terms of which we wish to 
expand the system hamiltonian, we can now rewrite our initial expression for the Bose 
fleld operator, Eq. (178), in terms of this new operator. Prom Eq. (183), we have* 

^(r, t) = ^/No [aj] 5^(r, t) . (185) 

The Bose fleld operator can now be written in the following number-conserving form 

*(r, t) = </.^(r, t)ao + ^/7Vo [aj] 5^(r, t) , (186) 

where we note that creation and annihilation operators do not commute, i.e. [oq]"^ 
oo, and {6'^{r,t)) = 0. We now substitute this into the system hamiltonian, Eq. 
(17), and consider contributions of an ascending number of 5^{r,t) operators. We 
also deflne the new effective interaction parameter ^jv via gN = gNo (rather than gN 
discussed in Sec. 5.2.1). To illustrate the procedure, we give here explicitly the flrst 
two terms, with calligraphic notation highlighting the fact that they now depend on 
the new fluctuation operators S^{r,t). So, 

Wo= y" dr[0^(r,i)]*M'^(r,t)aJao+| J dr |0^(r, i) f ajajaoao 
y dr[^^{r,t)r 

0"(r,t), (187) 





Here it becomes apparent how use of an operator under the square root (i.e. v A^o, or V N) would 
complicate matters considerably, as one would need to carefully consider the inverse of such operators, 
and corresponding operator ordering in the above expression. 
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where we have used the fact that ajaoao^o = a\{aQal — l)ao — {o^qClqY ~ (t*o'^o) = 
Nq - Nq. Similarly, we find 

Wi= 1 dr[,/.^(r,ir/ioaJ (^/7Vo[aJ]-l5^(r, i)) 

+ gj dr|</.^(r,t)|'[<^^(r,i)]*aJaJao(0Vo[aJ]-i^^(r,i)) 
+ h.c. 

= rfr[</,^(r,t)]*|Ao + 5iv |</)^(r, f ^^(r, i) 

+ h.c. , (188) 

upon re-expressing a\a\(XQ[al^^^ — a\{aQa\ — l)[a\]^^ — a\aQ — 1 = TVq — 1 . In Eqs. 
(187)-(188) we have divided all contributions involving a single condensate number 
operator, Nq, by its mean value {Nq) = Nq. The above analysis shows that Tip is 
of order A^O) and iii is of order ^/NQ. One can express all higher-order terms {ii.2, 
Tis and I-L4) in a similar manner, with each subsequent contribution Tim+i containing 
an additional factor of {1/\/Nq) and an additional fluctuation operator [5^{r,t)]^^^ 
compared to Hm (0 < rn < 3) - for detailed expressions see [161]. This justifies 
our earlier claim that we can expand the system hamiltonian systematically in terms 
of decreasing atom numbers, with the approach presented in this section effectively 
corresponding to an expansion in the ratio (iVjvc/-^o)- 

Following the methodology introduced in the number-conserving treatments of 
Sees. 2-4, we now proceed with suitable approximations for Hs and 7114, in order to 
obtain a more 'manageable' number- conserving hamiltonian. Although we cannot fully 
justify such approximations in the brief subsequent discussion, Gardiner and Morgan 
have argued in detail [161] that these approximations lead to the lowest non-trivial 
order which facilitates a consistent finite temperature description for a fixed finite 
number of atoms. 

Firstly, we entirely ignore the H4 contribution, and, for consistency (since the 
omitted contributions are of similar order of magnitude [78]) we simultaneously impose 
quadratic mean-field approximations - analogous to those of Eq. (53) - for the products 
of three fluctuation operators as 

[5^(r)]tj^(r')<5^(r") « ([(^^(r)]t5^(r'))<5^(r") 

+ ([,5^(r)]t,5^(r")>(5^(r') + [5^(r)]t(5^(r')^^(r")) (189) 

Simultaneously, we replace VVo by Nq in 'H2 and Hs. The 'dominant' contributions Hq 
and Hi should however be treated more accurately. In these expressions, we set 

NQ = NQ-Jdr {[5^(r)]t5^(r) - ([<5^(r)]tj^(r))} , (190) 

which physically corresponds to setting the number fluctuations between the 
condensate and the non-condensate to be equal and opposite [161] . 

The above approximations lead to a quadratic number-conserving hamiltonian. 
This generates the following set of non-local equations for the condensate and the 
non-condensate, which can be thought of as the number-conserving analogues of the 
HFB equations. The condensate energy evolves according to 

i/i^<^^(r,i) = [hQ{r,t) - A(i)] .^^(r,i) 
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+ 5(iVo(i) + A7Vo)|(/.^(r,t)P(/.^(r,t) 

+ 2gh{v, i)<^^(r, t) + 5m«(r, i) [<^^(r, t)] * - /(r, <)(191) 

The quantity f{r,t) appearing above is related to the dynamical interaction between 

the condensate and the non-condensate [165]; more specifically, it refers to the back- 
action of the (changes in the) normal and anomalous averages on the temporal 
evolution of the condensate mode, and is defined by 

f{r,t)= J dr'<?|0^(r',t)f {n(r,r',t)0^(r',t) 

+ m«(r,r',t)[</>^(r',t)]*} . (192) 

The physical justification for such a term is to ensure orthogonality is maintained 
throughout the system evolution. Moreover, AA^o = ((-^o) — {No)^)/No - 1 
is a typically small contribution from statistical fiuctuations included here for 
completeness. The number-conserving versions of the normal and anomalous 
averages are respectively given by h{r,r',t) = {{S^ {r,t))^6^ {r' ,t)) and m(r,r',t) = 
{d^{r,t)6'^{r',t)). The fluctuation operators evolve according to 

d ( 5^{v,t) 



C^{r,t) M^(r,t) \ / ^A^(r,t) 

M^{v,t)Y -[£^(r,t)l* J (<5^(r,t))t 



(193) 



where we have defined the corresponding generalized number-conserving expressions 
£^(r, t) = ho{r, t) - \{t) + gNo |0^(r, t) |' 

+ gNoQ{t)\^^{r,t)\'Q{t), (194) 
M''{r,t)^gNoQ{t) [0^(r, t)] ' 0*(t) . (195) 

One can make two immediate observations which are a direct consequence of 
orthogonality and number-conservation: (ij firstly, these expressions explicitly include 
the time-dependent projection operators Q(t) of Eq. (176), which guarantee that the 
orthogonality of the considered subspaces is maintained during the entire evolution, 
(ii) secondly, the derived equations of motion are highly non-local, with double 
projectors implying the following relation [166] 

Q(i)|</.^(r,t)|'Q(i)5^(r,i) = 

J y"dr'dr"Q(r,r',t)|</.^(r',t)|'Q(r',r",i)5^(r",i) . (196) 

Having obtained the full evolution equations for condensate and non-condensate 

operators, Eqs. (191) and (193), we can also investigate the frequencies of collective 
excitations in the presence of a dynamical thermal cloud, a regime in which this 
theory has been applied to with remarkable success (see Sec. 7.1). This is achieved 
by simplifying the full dynamical evolution in the context of a perturbative linear 
response analysis, analogous to that presented in Sec. 4.3.1. In particular, we consider 
perturbing the system by an external force which induces the mean fields to oscillate 
around their static values, i.e. 4>^{r, t) — (f)^ (r) + (5(/)^(r, t), n(r, t) ~ no(r) + (5n(r, t), 
and fh^{v,t) = rhQ{r) + 6-m{r,t). As before, we obtain the linearized equation of 
motion for the condensate fluctuations 6(f)^{r,t). In order to study the response of 
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the system to an external perturbation, we Fourier transform the temporal variables 
of (5(6^(r,t) and ((5(/)^(r, t))*, and expand the resulting expressions in terms of 
quasiparticle amplitudes via [166, 167, 168, 169] 



The response amplitudes bi (w) appearing above contain information on the condensate 
density fluctuations. These amplitudes are proportional to the response function 
which, in general, is made up of two contributions: a static one, arising from 
mean fields, and a dynamic one, associated either with direct excitation from the 
perturbation or indirect excitation from fluctuations induced in the other mean fields 
as a result of the perturbation. At low temperatures, the former contribution generally 
dominates, with the condensate excited directly from the applied perturbation. 
However, under certain conditions at finite temperatures, the perturbation may 
predominantly excite the non-condensate part, whose response to the perturbation 
leads to a secondary mechanism for the excitation of the condensate (see Sec. 7.1). 

5.4.. Brief Summary 

While mean field theories have been shown to produce generally very good results, their 
formulation is still based on the assumption that the condensate can be described as 
a macroscopic classical entity, thus crudely discarding its operator nature, which is 
generally 'allowed' for systems with a large number of atoms. Such treatment however 
automatically precludes the description of low- dimensional gases, where the simple 
picture of a coherent condensate is replaced by a so-called 'quasi-condensate', which 
can be thought of as a condensate which exhibits large phase fluctuations; an ab initio 
modified mean field theory has been formulated to deal with such cases at equilibrium. 

Moreover, the formal development of the preceeding sections where the system 
hamiltonian is separated into contributions according to the number of non- condensate 
operators appearing in them, can actually be generalized rather straightforwardly to 
the number- conserving case. To do this, one must first ensure the orthogonality be- 
tween condensate and excitations; this can be achieved by expanding the Bose field 
operator as ^{v,t) = (f>^{r,t)aQ + d{r,t), where ciq denotes an operator annihilating 
a condensate atom which is defined as the projection of the full Bose field operator 
^(r,t) onto the condensate state (j)^{r,t) (where N stands for Number-conserving) ; 
the non-condensate operator 6{r,t) is specified by the orthogonal projection. Identi- 
fication of the condensate wavefunction as the state with the largest eigenvalue via 
the Penrose-Onsager criterion, combined with the desire for the fluctuation operator 
to scale as the number of non-condensate atoms, leads to the definition of an ap- 
propriate non-condensate operator {r,t) = {aQ/\/N)S{r,t), which is guaranteed to 
satisfy {S^{r,t)) = 0. and in term,s of which the system hamiltonian can be expanded 
in a number-conserving fashion. Following pretty much the same procedure as in the 
mean-field treatments, the full system hamiltonian is now expanded in contributions 
of ascending orders of 5^{r,t); as each such term, additionally contains a prefactor of 
I/^/Nq, such an expansion can in fact be visualized as an expansion in powers of the 
ratio of non-condensate to condensate atoms; this justifies the perturbative treatment 
of hamiltonian contributions {H.^ + H.^) containing three or four number- conserving 
non- condensate operators S^{r,t). This approach leads to a finite temperature Gross- 
Pitaevskii and corresponding Bogoliubov equations which have similar structure as 




(197) 
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their mean field counterparts, hut also contain additional non-local terms arising from 
the orthogonality requirement. While this theory has been successfully implemented 
to the temperature dependence of the shift of elementary excitation frequencies in the 
linear response limit, a full self- consistent treatment of the thermal cloud is lacking at 
present. 

This concludes our perturbative treatment of generalized mean-field approaches 
and the related number-conserving formalism, on which this Tutorial has been largely 
based. The treatments presented thus far are summarized in Fig. 9, which makes a 
systematic classification of all previous approaches in terms of how the thermal cloud 
is treated, which hamiltonian contributions and (where appropriate) which generalized 
mean fields are included in each theory, what is the excitation spectrum predicted by 
each model, and what additional (e.g. many-body, low-dimensional) effects can be 
described by each such approach. 

The most advanced current approaches presented from each 'class' of models, 
namely the symmetry-breaking treatment of Sec. 4.4 and the number-conserving 
approach of Sec. 5.2.2 lead both to a successful numerical implementation and to a very 
good agreement with ultracold gas experiments at finite temperatures for a wide range 
of experimental conditions. However, it would be misleading to conclude this Tutorial 
without describing alternative rather distinct approaches that have been employed 
to study trapped Bose gases at finite temperatures. Such approaches, presented in 
Sec. 6, are actually required to accurately tackle various experimental issues, such as 
the onset of condensate growth from a purely thermal cloud, and fluctuations in low- 
dimensional systems. Most of the approaches presented below, which are based on 
very different theoretical formulations, are simply appropriate generalizations of the 
ideas of earlier sections, and can thus be shown to reduce to theories already discussed. 
For example, the successful approach of Zaremba, Nikuni and Griffin described in Sec. 
4.4.3 arises as a special case of the approaches of Stoof (Sec. 6.2.2) and Gardiner-ZoUer 
and co-workers (Sec. 6.2.3) in the appropriate limits. It is important to stress that 
while such approaches are presented last here, as appropriate from a pedagogical point 
of view, they have been actually derived over a large number of years, with some of 
the most important results preceeding those of alternative treatments of Sees. 4.4 and 
5.2.2. 

Each of the topics mentioned below would merit a Tutorial in its own right, 

and, for example, a Review on Classical Field Theories has been published recently 
[53] . For the purposes of this Tutorial, the subsequent treatment is limited to a brief 
overview, highlighting the main physical ideas, final equations and applicability of 
such approaches. Numerous references are given to aid the reader who may wish to 
pursue further studies on these topics. An excellent set of presentation transparencies 
(with additional references) on these topics can be found in the websites of two 
recent meetings on finite temperature models for Bose-Einstein condensation held in 
Sandbjerg (Denmark) [170] and Heidelberg (Germany) [171] in the summer of 2007. 

6. Alternative Beyond Meein Field Approaches 

The approaches presented below are explicitly number-conserving, and are already 
widespread in the study of ultracold atomic gases (see also [52] for a related review of 
such approaches). 
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6.1. Classical Field Approaches & The Projected Gross-Pitaevskii Equation 

In our preceeding discussion we have argued against the validity of the Gross-Pitaevskii 

equation at finite temperatures, on the basis that it ignores both mean field effects 
and dynamical couphng to the thermal cloud. However, the Gross-Pitaevskii equation 
is actually an equation for a classical field; as such it should be able to describe all 
classic;al aspec;ts of a finite temperature system of ultracold gases. For example, it is 
well-known from the classical theory of electromagnetic radiation that the Rayleigh- 
Jeans model provides a good approximation to the full quantum system for modes 
with energies less than kBT. provided all suc;h modes are 'highly occnipiecr, such 
that the Bose-Einstein distribution function of Eq. (71) can be well approximated by 
the 'classical' expression /(e^) « {kBT)/ei. Formalisms based on this approximation 
are termed 'classical field' approaches and have recently been reviewed in [53]. The 
discussion of this approximation in the context of ultracold Bose gases has been driven 
by work of Svistunov, Kagan and Shlyapnikov on the formation and dynamics of 
Bose-Einstein condensation [172, 173, 174, 175]. Such an approximation was also 
used early on in other contexts, such as the electroweak phase transition [177] and 
the equilibration of a Bose gas to a superfluid state [178], with a qualitative two- 
dimensional study of evaporative cooling performed in [179]. 

To study an ultracold Bose gas, one seeds the initial state of the system 
with arbitrary initial conditions. Typically one expands the initial wavefunction 
in appropriate eigenstates ipk(j) via </>(r,t = 0) = Cfe(^/j(r); in doing so, the 
amplitudes and phases are appropriately chosen under the constraint of fixed total 
atom number and energy: for example, for the homogeneous gas, one typically chooses 
(^fe(r) = e*'^ '', with the populations |cfe|^ chosen such that the distribution is as fiat as 
possible, and the phases of Ck are selected randomly [175, 180]. Such an initial state 
is then propagated by the GPE of Eq. (29). Due to its intrinsic nonlinearity, this 
equation mixes different modes and relaxes rapidly to a classical thermal distribution. 
In fact, the precise initial conditions are largely irrelevant, as they are lost after a 
short temporal evolution, and it is therefore not even important to set the initial 
conditions to be equilibrium ones for the particular system. In these simulations, 
the temperature of the system is not set a priori, but it is actually subsequently 
determined by extracting the temperature upon fitting the number occupation of 
the relaxed system with a classical thermal distribution. The parameter (j){v,t) in 
such simulations does not simply model the condensate, but actually it describes the 
entire multi-mode 'classical' atomic gas, thus enabling the study of various parameters, 
such as correlation functions. Clearly such simulations require an upper energy 
(and momentum) cut-off to ensure that the system remains in the 'classical regime' 
throughout the simulations [53]. 

Although the use of a numerical grid in performing simulations imposes itself 
a cut-off, the most accurate implementation of this approach explicitly introduces 
a projector into the GPE, to ensure that only momentum-conserving processes 
consistent with the large occupation approximation are considered, and no numerical 
'aliasing' arises [181, 182]. The projector used in this theory is defined by P{</)(r,t)} = 
"Y^keC^kiy) I '^''Vfc(i'')?^(i'') where k is restricted within the coherent (classical) 
region C [182]. This projector is only required in the nonlinear term, and leads 
to the replacement of the nonlinear term \4>\'^4' appearing in the original GPE by 
P{\(j)\^(l)}. The resulting equation, introduced by Davis, Morgan and Burnett is 
known as the Projected Gross-Pitaevskii Equation (PGPE) [180, 182, 183]. While 
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N /N = 0.93 ■ N /N = 0.45 ■ N /N = 0.02 
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Figure 10. (colour online) Typical thermalized (equilibrium) images of a 
classical field consisting of a fixed number of 2, 000 atoms at three different 
energies (i.e. temperatures). Plotted are the density profiles |</>(a;, y, 0)p of an 
anisotropic (ajj^/ojz = \/8) trapped 3D Bose gas arising from a single run of 
the PGPE, with colour representing the atomic density (plotted on logarithmic 
scale - black/blue: zero density, red/dark brown: maximum density). The 
enhancement of fluctuations at higher temperatures corresponding to a lower 
condensate fraction Nq/N (indicated on figure) is evident. (Images provided by 
Matt Davis - see also [f83].) 



the implementation of the projector is very straightforward in the homogeneous gas 
[180], its extension to the trapped case has only been discussed recently [183], by 
expanding the classical fields in harmonic oscillator eigenstates. Thermal fluctuations 
are included into the treatment, as evident from characteristic single-shot images 
of the density profiles of a trapped gas at three different temperatures, which is 
plotted in Fig. 10. However, other results such as condensate fraction and correlation 
functions which require suitable averaging over the fiuctuations, are evaluated by 
replacing ensemble averages by time-averages, based on an ergodic hypotehsis (see 
e.g., [53, 183]). Although such theories appear to work well in diverse contexts, an 
extension of the PGPE has also been formulated, in which the coherent region is 
explicitly coupled to a heat bath of non- condensate atoms [119]; in this treatment, 
additional terms are introduced into the PGPE to describe the coupling of the modes 
in the coherent region to a heat bath of non-condensate atoms. This idea of splitting 
the description of the system into low- (coherent) and high-lying modes (above a 
carefully selected cut-off) essentially constitutes a (multi-mode) generalization of the 
ZNG formalism, and will be further discussed in Sees. 6.2.2-6.2.3. 

There have been numerous applications of the classical field method [175, 176, 184, 
185, 186, 187, 188, 189, 190, 191] and of the PGPE [180, 182, 183, 192]. The PGPE 
is a non-perturbative method and has also been used to study the process of non- 
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equilibrium condensation [176, 193], to determine the shift in the critical temperature 
Tc [192], and to study spontaneous vortcx-antivortcx pair production in quasi-2D gases 
[194]. Despite their success, classical field approaches discard quantum fluctuations 
hy construction, and therefore more accurate treatments may be needed in certain 
contexts - as developed below. 

6.2. Stochastic Methods 

We thus discuss next various schemes which introduce quantum fluctuations into the 
theoretical description in different levels of approximations. 

6.2.1. The Truncated Wigner Approximation: Closely related to the above approach 
is the so-called Truncated Wigner approximation [195]. In this approach, quantum 
fluctuations are approximately included in the formalism as follows: in particular, the 
system evolution is monitored as before by the usual (Projected) Gross-Pitaevskii 
equation, only the initial conditions have been modified to additionally contain 'the 
right amount' of quantum noise, i.e. quantum fluctuations of half a particle per mode 
(on average), such that they appropriately sample the full Wigner function. The effect 
of quantum fluctuations is thus introduced here by averaging over numerous distinct 
numerical realizations. 

This theory can be put on firm ground as follows (an excellent detailed discussion 
of this can be found in [181]): Consider a system described by a density operator 
p{t). To investigate dynamical effects, instead of studying the evolution of the full 
density operator, one can equivalently monitor the evolution of a suitably-constructed 
quantum quasi-probability distribution function [154, 155, 156], known as the Wigner 
function, W. This is constructed in terms of a suitable integral depending on p{t) and 
on the eigenvalue, a, of the single-paxticle annihilation operator in a coherent state 
\a) (in the general case of a multi-mode field applicable here, this would be replaced 
by multiple integrals). One must thus consider the equation for the evolution of the 
Wigner function, which acquires a relatively straightforward form of a nonlinear partial 
differential equation in time. Since the phase space occupied by the Wigner function 
is quite large, one would instead prefer to track the trajectories of single realizations of 
the system through phase space, which would be described by some sort of stochastic 
differential equations. Knowledge of a large number of such trajectories, would then 
essentially enable the reconstruction of the Wigner function, and thus provide details 
of the system evolution. 

Consider for the sake of our current discussion a single-mode system with 
a corresponding Wigner function W{a,a*,t) (generalizes trivially to multi-mode 
systems): the full evolution of this Wigner function is known to contain, among 
other terms, third-order derivative terms of the form dW/dt oc [d^ /d'^ada*]{aW) — 
[d"^ /dad'^a*]{a*W), which do not facilitate an immediate correspondence of this 
equation to a stochastic differential equation [154, 155]. Fortunately, it turns out that 
these terms are quite small provided essentially that all modes of the system are highly 
occupied (/, ^ 1) (for details and precise conditions see [196]). In this limit, one can 
'truncate' the exact equation of motion by discarding such terms. The evolution of 
the Wigner function is thus mapped onto the evolution of a field, whose dynamics 
is governed by an equation formally identical to the Gross-Pitaevskii equation of 
Eq. (29) (and the same arguments regarding the projector apply as in last section) 
[181]. Having made this identification, one can now routinely solve the (P)GPE as 
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usual, and simply seed the initial states with noise, such that spontaneous processes 
also become feasible [181]. To further justify the above, we note that Polkovnikov 
recently performed a systematic perturbation theory in quantum fluctuations around 
the classical system evolution, obtaining the GPE to lowest order, and the truncated 
Wigner approximation to next order [197]. 

Although the truncated Wigner method has been used extensively in 
quantum optics [198], the first numerical implementation of the Truncated Wigner 
approximation in utlracold gases was performed in 1998 in [195]; a sampling 
technique based on the number-conserving Bogoliubov theory [157, 158] was then 
subsequently analysed in [196, 199] (see also [200]). The Truncated Wigner has 
since been used to study various phenomena, including colliding condensates [201], 
condensate reflection from a steep barrier [202], three-body recombination processes 
[203], collapsing condensates [204], atom interferometers [206] and optical lattices 
[207, 208, 209, 210, 211]. 

Both treatments presented above describe an ultracold atomic gas by a single 
equation for the highly-occupied modes of the system. At low temperatures, high 
mode occupation is restricted to the lower part of the energy spectrum, corresponding 
to the condensate and the states dressed by their proximity to the condensate. 
Thermalization is actually achieved by the nonlinear mixing of different modes. Due 
to their similarities, we should now briefly comment on the physical interpretation and 
validity of the above approaches. In the classical field approaches, one assumes that 
all relevant modes of the system are highly-occupied, such that thermal fluctuations 
entirely dominate over quantum fluctuations. This automatically restricts its validity 
above a minimum temperature (which is nonetheless relatively small), with the same 
method clearly also applicable above the critical region, where the gas is classical 
anyway. On the other hand, the Truncated Wigner method only differs from the 
PGPE in that it additionally includes numerical noise in the initial conditions of the 
simulations, which is aimed at mimicking quantum fluctuations. It would thus be 
natural to assume that the validity regime of the Truncated Wigner is the same as 
that of the PGPE, with an additional extension to the very low temperature regime - 
this is however not the case. 

As mentioned earlier, the Truncated Wigner method is based on an approximate 
equation for the evolution of the full Wigner quasi-probability distribution (which 
discards certain contributions), an approximation which should generally be valid as 
long as all relevant modes of the system are highly occupied. Importantly, in the 
Truncated Wigner method, noise is only included in the initial conditions in a manner 
which on average corresponds to the addition of half a particle per mode. However, 
evolution of this initial quantum distribution by the classical (P)GPE results in the 
initial quantum distribution thermalizing into a classical fleld at a slightly higher 
temperature than the initial distribution, thus giving rise to unphysical heating [196] . 
The net result of this is to limit the applicability of the theory either to relatively short 
times (before any such classical thermalization takes place), or to the low-temperature 
regime defined by the condition that the maximum energy of the Bogoliubov modes 
should not exceed a few times the thermal energy fc^T [196]. As a result, the Truncated 
Wigner method is more useful for studying situations where the quantum processes 
largely dominate over thermal effects. On the other hand, one could interpret the 
(P)GPE as the thermodynamic limit of the equation for the dynamics of the Wigner 
function, for which spontaneously initiated (as opposed to stimulated) processes are 
negligible. In other words, the Truncated Wigner method yields a more accurate 
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description of relatively low temperatures, where quantum effects provide a noticeable 
contribution, with higher temperature generally described quite well by the classical 
field method, or the PGPE. 

How can such an unphysical thermalization arising within the Truncated Wigner 
method be avoided? For this purpose, it is important to note that the low-lying modes 
we have been discussing so far are actually coupled to a range of higher-lying 'thermal' 
modes which are occupied more sparsely, and hence above the chosen cut-off. In fact, 
coupling the previously-considered part of the system to a thermal gas provides the 
necessary irreversibility for the system to relax, while simultaneously guaranteeing 
it relaxes to the correct equilibrium. This additional coupling gives rise both to 
dissipative terms and to dynamical noise terms which modify the Gross-Pitaevskii 
equation to an appropriate Langevin, or Stochastic Gross-Pitaevskii equation. 

Next, we describe two seemingly very distinct approaches, formulated respectively 
by Stoof [127, 128, 212], and Gardiner, ZoUer and co-workers [130, 213, 214, 215, 216], 
which are however essentially both based on a physical description of the system in 
terms of an appropriate probability distribution function. Such treatments combine 
the ideas of the formalism of Zaremba, Nikuni and Griffin (Sec. 4.4.3) whereby the 
gas is split into two dynamically coupled subsystems, with the inclusion of quantum 
fluctuations discussed previously. This is achieved by a generalization of the earlier 
picture of splitting the system into condensate and non-condensate contributions, to 
a separation between low- and high-lying modes. The low-lying modes forming the 
coherent region refer to the condensate and modes affected by it due to thermal and 
quantum fluctuations and are treated exactly. Such modes are dynamically coupled to 
the incoherent part of the system (high-lying modes) , as in the spirit of the approach 
discussed by Zaremba, Nikuni and Griffin, only their evolution is now stochastic, i.e. 
it includes explicit noise terms to account for fluctuations [128, 215, 216, 220]. This 
description arises by mapping the effective equation for the probability distribution 
of such low-lying modes onto a stochastic equation for their evolution. Although 
the non-condensate modes are usually treated as a thermal particle reservoir, the 
full description models the high-lying modes by an appropriate quantum Boltzmann 
equation. This set of equations, which in the authors' opinion constitutes the 'next 
generation' kinetic theory of ultracold bosonic gases, is mathematically closely related 
to the treatment of Sec. 4.4.3, although the interpretation of (/'(r,f) is now distinct, 
as it constitutes the 'order parameter' of the (multi-mode) system, rather than the 
condensate wavefunction. These theories can be shown to reduce in the appropriate 
limits to either of the self-consistent Gross-Pitaevskii-Boltzmann, or 'ZNG' theory 
(and thus to all more elementary mean field theories), the classical field method and 
the truncated Wigner approximation. 

6.2.2. Stoof s Non- equilibrium Theory: Stoof derived an equation for the evolution 

of a Wigner probability distribution P[(/)*.0;t] based on functional integration 
techniques [127, 128, 129, 212]. This distribution expresses the probability of the 
system to be in a coherent state \(j){r);t), derived from the vacuum state |0) by 
\(f>{r);t) = exp{/ dr(t>{r)^^ {r,t)}\{)) . By expanding the density matrix in coherent 
states, the probability distribution is cast in a form which requires the calculation of 
a functional integral containing the quantity \{(j);t\(j>o',to)\'^ = {(j); t\(j)o; to) {(j>o; to\(j); t) . 
Using standard techniques, each of these matrix elements (0; t|c/)o; to) can be written 
as a ^path integral' [57] over all complex field evolutions tp{r,t) (with ip{r,to) = ^o(r) 
and ^/'*(r, t) = 4>*{r)), via / d[ip*]d[tp]exp {iSEFF[''P* ,ip]/h}, where S'eff is the effective 
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action of the system (essentially analogous to the choice of an effective hamiltonian 
in our preceeding forinalisin). The determination of | (0; f|(;6o; /:o) | requires one to 
study all possible paths that the evolution of the field may follow from some time t 
back to the initial time <o (evolution backward in time), and then back again to t 
(forward evolution) , leading directly to the introduction of the Keldysh contour [46] , 
over which the paths are to be calculated. Such an approach is explicitly number- 
conserving, and typically discards any details of the initial quantum state (Markov 
approximation). After performing a standard transformation in the field variables 
ij} to explicitly separate the semi-classical dynamics from the effect of fluctuations, 
this treatment leads to a quadratic effective action in the fluctuations, in which the 
interactions have been renormalized to include many-body effects. The evolution 
of the probability distribution takes the form of a Fokker-Planck equation (see Eq. 
(198)), which is the generalized analogue of the 'quantum' evolution described by 
the truncated Wigner approximation. Fluctuations around can be systematically 
calculated, with correlation functions of any order obtained from suitable moments 
of P[(j)* ,<j); t]. An excellent non-technical discussion of this approach can be found in 
[218], whereas a more mathematical paedagogical presentation is given in [219]. 

The system may be separated into two 'components' by expressing the probability 
distribution of the system as P [0*, (j); t] = Pq [$*,$; t] Pi [(j)'* , (j)' ; t], which amounts to 
a Hartree-Fock-type approximation. This procedure leads to two coupled equations: 
firstly, substituting this ansatz into the full Fokker-Planck equation and integrating 
over the non-condensate degrees of freedom (j)' leads to a Fokker-Planck equation (given 
by Eq. (198)) for the low-lying modes of the system, which constitute the 'coherent 
region'. Performing the related integration over the condensate contribution $ leads 
to the corresponding Fokker-Planck equation for the non-condensate, whose 'semi- 
classical' treatment yields a Quantum Boltzmann equation (Eq. 203)) describing the 
incoherent system dynamics. 

The low-lying modes of the system are described by the Fokker-Planck equation 
[128, 129] 

i^l^Po [$*,$;*] = 

"/ ^""^ [h^-IJ.{t)-iR{v,t)+g\^{vf] $(r)Po[**,*;t] 

'K - m + iR{r, t)+g mr)f] $*(r)Po [$*, t] 



dr- 



(5$*(r) 

'^ ^^(r)l-(r) ^^(^^^^"(^'^^^°[^*'^'^J ' 

where 5/5f^^*\r) denote functional derivatives [69, 57]. Let us now briefly interpret 
this equation: Both iR{v,t) and S^(r,t) arise from the coupling to the higher-lying 
modes (the 'reservoir'). Let us flrstly ignore the last contribution of Eq. (198). If 
we were to additionally set iR — 0, then the above equation would be equivalent to 
the GPE for ^>(r, i) (flrst line), and the corresponding evolution for $*(r,t) (second 
line), which would thus be describing the classical evolution of the low- lying modes 
of the system. The iR, term is precisely the contribution that was included in our 
perturbative beyond-HFB approach and describes the transfer of particles between 
the coherent part of the system and the higher-lying (thermal) modes of the system. 
It is given by Eq. (153), upon a slight redefinition whereby we drop the Vc term from 
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the momentum-conserving delta functions and absorb the kinetic energy into /Xc, such 
that Ec of Eq. (149) is replaced here by = /i^; moreover the non-condensate energies 
ei of Eq. (148) are replaced by ei(r,p) = (|pi| V2rn + Kxt(ri) + 2g{\^{vi)\'^)). The 
explicit inclusion of iR modifies the 'reactive' GPE to the dissipative GPE of Eq. 
(152). The latter contribution in Eq. (198) contains the so-called Keldysh Self-energy, 
S'f^(r,t), given by [128, 129] 

v^r +\- _ Y ^\ 2 / <ip2 f dp3 f dpi 
1. [r,t)- ^[j,)9 J (2^^)3 J ^2nnr J {2nhr 

X {2'Ktl)^5 (P2 - P3 - P4) S {Sc + £2 - £3 - £4) 

X [{N2 + 1)N3N4 + N2{N3 + 1){N4 + 1)] , (199) 

where Ni = N{ei) denote the thermal populations. corresponds to fluctuations 
associated with the same collisional processes described by iR (i.e. incoherent collisions 
between condensate and non-condensate atoms), so the last contribution of Eq. 
(198) should also be included whenever collisions are included into the treatment 
via iR ^ 0. Its inclusion maps Eq. (198) onto a stochastic partial differential 
differential equation for this equation is essentially equivalent to the dissipative 
Gross-Pitaevskii equation discussed in our second order mean field approaches of Sec. 
4.4 (see, e.g. 'ZNG' Eq. (152)) but with the further addition of an appropriate noise 
term. 

However, both the dissipation iR{r,t) and the self-energy S^(r,i) depend 
implicitly on $(r,t) through their dependence on the condensate energy £0 = ^,0 + 
g\^{r,t)\'^ and through the induced change in the non-condensate factors N^. Thus, 
the corresponding stochastic equation is a Langevin equation with multiplicative noise 
and a prefactor with a complicated dependence on $(r,i) [129], which is hard to 
simulate numerically [129]. To reduce this equation to a more manageable form, we 
restrict our discussion to the regime near equilibrium, for which one can show that 
the Fokker-Planck equation reduces to the following stochastic equation [129, 220] 



l+^;iE^(r,t) 



(ho + gmr,t)\^-n)^r,t) 



+ V{r,t) (200) 

with the noise term having gaussian correlations of the form 

(ry*(r,t)r/(r',t')> = «(fiV2)S^(r, t)(5(f - f')'5(r - r') . (201) 

How does such an equation come about? From the formulation of the theory, close 
to equilibrium where the thermal cloud can be described by the Bose distribution 
function, one can derive a fluctuation-dissipation theorem linking the strength of 
the fluctuations, HT,^ , with the dissipation, iR. Both of there depend on the 
energy Sc for promoting a condensate atom (located at position r) out of the low- 
lying part of the system, with their mathematical relation given by iR{r,ec) = 
— {h/2)T,^ {r,ec)[l + 2A^(£c)]~^. Using now the classical field approximation N{ec) ~ 
{UbT) I (£c ~ ii) valid for the low-lying modes of the system which are highly populated 
(iV(£c) > 1), one can re-expn-ss [1 + 2N{£,)]-^ « (/^/2)(£c - m)- Finally we note 
that, within the present formalism, Sc is actually an operator in the configuration 
space of the order parameter in contrast to the semi-classical treatment of Sec. 
4.4.3; this leads directly to the approximate expression for the dissipative term 
iR{r,t) « — (/3/4)/iS^(r, t)[^o + 9\^{''^,t)\^ — A*], which has already been implicitly 
assumed in Eq. (200). A plot of the position and temperature dependence of the 
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Position 



Figure 11. (colour online) Typical position and temperature dependence of the 
self-energy for a trapped one-dimensional Bose gas {u)^ = 2tt X lOHz, *^Rb). 
Plotted is the scaled self-energy filKE {z)\ / A appearing in Eq. (200) (where 
/3 = 1/kgT) against position (measured in harmonic oscillator units) for two 
different temperatures T = 400nK (solid black, top) and T = lOOnk (dashed red, 
bottom). 



self-energy , E^(z), for a harmonically trapped one-dimensional Bose gas is shown in 
Fig. 11. 

Having already identified the origin of the different contributions, we can now 
re-interpret the quantity ^(r,t) appearing in Eq. (200) as the order parameter of 
the system, since it describes not only the condensate, but also incorporates thermal 
and (in its full version) quantum fluctuations, i.e. accurately describes the low-lying 
modes of the system [128, 220]. We can now identify the important new physical 
feature introduced by this equation. All mean field approaches discussed earlier had 
the fundamental restriction of no spontaneous initiation, i.e if the condensate mean 
field were initially zero, it would remain so at all subsequent times. Contrary to that 
picture, the presence of the noise term in Eq. (200) provides a 'quantum mechanical 
seed' to initiate condensate growth. (In fact, the critical region has been studied in 
more detail by Stoof, who showed explicitly that the system acquires the necessary 
irreversibility required for the onset of condensation [127, 128].) 

The derivation sketched above relies on the quadratic nature of the effective 
action, and some comments are necessary here: Although the cubic and quartic 
contributions of fluctuations around (j) do not explicitly appear in the treatment, their 
effect has actually been taken into consideration in generating both the effective many- 
body T-matrix interaction and the coupling to the high-lying modes. This argument is 
essentially equivalent to the elimination of the anomalous average to yield many-body 
effects (Sec. 3.4), and the interpretation of the triplet correlations (S^SS) as providing 
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Figure 12. (colour online) Typical evolution predicted by the solution of the 
stochastic Gross-Pitaevskii equation in two dimensions. Top/Middle; Single-run 
results showing the density/phase, clearly highlighting the spontaneous vortex 
formation in the ab initio growth of a quasi-condensate coupled to a static heat 
bath. The quasi-condensate is formed in a harmonic trap with u)^ = = 
2-n X 200Hz, and contains approximately 65000 -^^Na atoms at equilibrium for a 
heat bath at T = 500nK. Bottom: Corresponding profiles obtained after averaging 
over a small number of runs (si 25), as might be done experimentally; in the 
latter profiles spontaneous singularities are largely washed out. (Images provided 
by Stuart Cockburn - see also [221]; note that in the above images the colourbar 
is redefined from column to column for optimal visualization.) 



the required coupling to the non-condensate (Sec. 4.4). 

To complete our treatment, we must also discuss the evolution of the occupation 
numbers N . This can be done by introducing the Wigner distribution via [128] 

j dr'e-P "-'/" (0' [v + r (r - ) = ^(P, i) + ^ , (202) 

and using the corresponding Fokker-Planck equation for the non-condensate 0' 
(obtained by integrating out the condensate degrees of freedom) to extract the 
evolution of {(f)'(j)'*). Comparison to the non-condensate populations /(p,r,t) 
introduced in Eqs. (143)-(144) reveals that the field (f>'{r,t) now contains both thermal 
and quantum fluctuations, with quantum-mechanical fluctuations of half a particle 
per mode added to the thermal occupation N{p,r,t) of each mode. Although such a 
noise contribution was previously introduced in the initial conditions in the Truncated 
Wigner approach, here the noise is dynamical, i.e. it is added onto <I>(r,t) at each 
time-step in the simulations [220]. The non-condensate populations evolve according 
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to a quantum Boltzmann equation. Making the usual gradient expansion (assuming 
the non-condcnsatc varies on a much larger lengthscale than the external trapping 
potential), one obtains [128] 

dN 

— + (Vp£) • (VrTV) - (Vr£) • (VpTV) = C^2[N] + CasfAT] . (203) 

The two coUisional integrals Ci2[A'^] and C22[^] appearing above arc practically 
identical to those discussed in Sec. 4.4.3 in the context of the 'self-consistent Gross- 
Pitaevskii-Boltzmann', or 'ZNG' theory. (Eqs. (157)-(159)). One important difference 
here is that they arc expressed in terms of the quantity $(r, f), instead of the 
condensate density 4'{v,t): 

X {2nh)^S (p2 - P3 - P4) <^ (£c + £2 - £3 - £4) 

X {2nHf [5{p - P2) - S{p - Pa) - S{p - p^)] 

X [{N2 + l)NsNi - N2{Ns + l)(iV4 + 1)] , (204) 

47r 2 f dp2 f dp3 f dp4 



r \m = — o2 f dp2 f dps f 

^ J {2T:hf J (27rfi)3 J (27rfi)3 



X {2iThf5 (p + P2 - P3 - P4) 5 (e + £2 - £3 - £4) 

X [{N + 1){N2 + l)N3N4 - NN2{N3 + 1){N4 + 1)] . (205) 

The stochastic approach presented in this section amounts to solving self- 

consistcntly the stochastic Gross-Pitacvskii equation, Eq. (200), with dissipation and 
noise respectively defined by Eqs. (153), (199) and (201), coupled to the dynamical 
Quantum Boltzmann Equation of Eq. (203) with the coUisional integrals of Eq. (204)- 
(205). The non-condensate density is obtained via n^c = J dp/{2TTh)^N{p,r,t) (as 
in Sec. 4.4.3). Applications of this approach (in suitable limits) can be found in 
[129, 151, 220, 222, 223, 224, 225] (see also Sees. 7.1-7.2). 

How arc the solutions of the Stochastic GPE to be interpreted? As the formalism 
explicitly includes dynamical noise, the main idea behind it is to extract information 
from an averaging over a large number of numerical 'runs' (ranging from a few tens 
to a few hundred - depending on required accuracy and regime of simulation). Due 
to the random nature of fiuctuations included, averaging over many runs is equivalent 
to averaging over different experimental realizations; such an approach yields, for 
example, smooth density profiles and correlation functions. It should however be 
remarked that as the Stochastic GPE describes the whole matter-wave field, it makes 
no distinction between condensate and thermal cloud, a common 'artificial' feature 
of all treatments presented in Sees. 2-5: To highlight this point the images of the 
'modified low-dimensional theory' presented in Sec. 5.1 were compared to results of 
the stochastic code for the same parameters (see earlier Fig. (8)); while the overall 
agreement was very good, the stochastic images automatically produce the total 
atomic density, rather than needing to sum over the two independently-determined 
(but coupled) sub-contributions. Thus simulations of the stochastic GPE can be seen 
as corresponding precisely to 'numerical experiments', and therefore an (artificial) 
partitioning of the atomic gas to a 'condensate' and a 'non-condensate' contribution 
should be performed in the same way as when analyzing experimental data - e.g. via 
bimodal density fits [61] (or via a manipulation of system correlations [224]). 
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At the same time, like with the Projected GPE discussed earUer, one can also 
extract useful information from a single run of the Stochastic GPE, which would 
resemble results from a single experimental realization - such analysis is relevant when 
looking at processes which are sensitive to individual runs, such as the collapse of a 
BEG [3, 129], or the visualization of images of phase and density fluctuations. For 
example, a single run of the Stochastic GPE in two dimensions (Fig. 12, top images) 
leads to the spontaneous appearance of vortices as the (quasi)condensate grows from 
the thermal cloud, with the location (and the number) of the vortices differing from 
realization to realization (see also [221]) - the spontaneous appearance of vortices 
in the context of a second order phase transition induced by a rapid quench was 
originally proposed in the cosmological constant [227, 228], and subsequently discussed 
in the context of an atomic BEG in [229]. However, as soon as one averages over a 
large number of independent realizations (runs) spontaneous excitations are rapidly 
'washed out', leading to much smoother density proflles (Fig. 12, bottom images), as 
also relevant for calculations of correlation functions [224]. 

A simplified form of the Stochastic GPE of Eq. (200) with time-independent 
occupation numbers in the non-condensate, and thus a time-independent self-energy 
E^(r) (see subsequent Eq. (216)) was used to study reversible condensate formation 
when cycling through the phase transition [226], thus providing the very first 
application [220] of such a stochastic equation to the study of ultracold gases. One of us 
(NPP) has been heavily involved in subsequent simulations, investigating fluctuations 
of one-dimensional Bose gases [151, 224], the growth of coherence of an atom laser 
[225], and quasi-condensate growth on an atom chip [223]. It is also important to 
remark that instead of solving the full SGPE numerically, one can resort to variational 
calculations. Such a technique was used to discuss collisional frequencies and damping 
rates of collective excitations, growth-collapse cycles in attractive condensates [129], 
and flnite temperature dynamics of a single vortex [222] . 

Stochastic Hydrodynamics: Gonsideration of the stochastic effects leads to 
a modiflcation in the corresponding finite temperature hydrodynamic equations 
discussed in Sec. 4.4.4. To study these, we use the Madelung transformation 
$(r, t) = i/n(r, f)e'^('''*^ directly within the effective action, which is thus re-expressed 
as S'eff['^, ^] in terms of the density n{r,t) and the phase 0{r,t) of the system. This 
procedure leads to the following two coupled stochastic equations of motion [129]: 



where the above noise sources (^„, have gaussian correlations of the form 



Here flc = Ho + 9n+ {l/2)mv'^ = -(fi^V^-yn)/(2mv^) -|- Vtrap + gn+ (l/2)mvf is 
the chemical potential of the condensate (the condensate kinetic energy contribution 
has been absorbed into the chemical potential here, contrary to the corresponding 
definition of Eq. (139), and /x is the chemical potential of the thermal cloud. 
Gomparing these equations to the hydrodynamic description given earlier, we find 
the following modifications [129]: (i) The 'continuity equation' acquires an additional 
drift term due to (ii) The equation for the condensate phase (from which the Euler- 
like equation for the superfluid velocity can be obtained via v(r,f) = {h/ni)V9{r,t)) 




(206) 



(207) 



{U^Mn{v',t')) = {UrMe{r',t'))^t-E^{r,t)Sir-r')6{t-t').{208) 
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acquires both a spatial diffusion-like contribution proportional to iT,^{v,t) due to 
collisions between condensate and thermal atoms, and a noise contribution inversely 
proportional to the square root of the density. 

Next, we briefly review the results of Gardiner, Zoller and co-workers who 
used a different methodology, which is nonetheless similar in spirit to the above 
treatment, and also yields a similar Stochastic Gross-Pitaevskii equation. Before 
doing so, we should however also highlight here the related work of Ramos and 
co-workers [230, 231, 232] on the non-equilibrium field-theoretic formulation and 
numerical simulation of the equilibration of a homogeneous Bose gas following a rapid 
temperature quench. 

6.2.3. The Gardiner- Zoller Kinetic Theory: The theoretical description of Gardiner, 
Zoller and coworkers [130, 213, 214] is also explicitly number-conserving and is based 
on techniques established in the quantum optics community, and described in various 
textbooks [154, 155, 156]. In their treatment, the system is also split into two parts: 
the first one, termed the 'condensate band' {Rc) corresponds to the low-lying modes 
of the system, which include the condensate and those modes which are affected by 
its presence; states lying above a particular energy, belong to the 'non-condensate 
band' {Rnc)- Note that, as before, the non-condensate band includes all modes 
wliic;li may be occupied during a collision, and should not be confused with the yet 
higher-lying modes which have been implicitly eliminated in order to introduce the 
effective interaction in the binary collisions. These two bands, apart from having their 
own internal dynamics, are also allowed to interact and exchange both particles and 
energy, pretty much as in the schematic of Fig. 5. The system hamiltonian is thus 
split into three parts, describing the contributions within each band and HNc)t 
and their respective interaction. Hint- One can write down an equation of motion 
for the density operator of the system, p, via [130, 213, 214] 



By making the assumption that there are many more atoms in Rnc compared 
to those in Rc, one may initially assume that the non-condensate remains practically 
unaffected by the interaction process, and can thus be treated as a 'reservoir'. In order 
to obtain the evolution of populations in the condensate band, the standard procedure 
then is to eliminate the reservoir by tracing over its degrees of freedom, so that one 
ends up with a density operator for the condensate, pc = TrNc(/5)- This is similar 
in spirit to the elimination of the non-condensate modes 0' performed by Stoof in 
order to obtain a Fokker-Planck equation for the low- lying modes of the system. One 
then uses the Laplace transform method (and also ignores interactions in the kernel 
and makes the Markov approximation [130]) to obtain the so-called 'master equation' 
for the evolution of p. This is a fully quantum-mechanical equation describing the 
evolution of the low-lying modes in contact with the non-condensate which is treated 
as a heat bath. The derivation of the resulting master equation is quite lengthy and 
will not be given here, with the final form of this equation presented in Appendix C.l. 
Suitable limiting cases of this master equation have already been used in numerical 
simulations [120, 233, 234, 235, 236, 237, 238, 239, 240, 241], and are briefly presented 
below. We should also mention here the related work of Anglin [242] who formulated 
a master equation for a single trapped mode of an ultracold weakly-interacting Bose 
gas. 




(209) 
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The 'Quantum- optical Master Equation': In its simplest application, the master 
equation can be used to generate rate equations for the average populations of each 
mode - as done in textbook discussions of laser theory [154, 155, 156]. The simplest 
approximation that can be used to study condensate dynamics, is to 'shrink down' 
the entire condensate band to simply one mode describing the condensate (as in 
the treatment of Sec. 4.4.3), and study its interaction with the reservoir of thermal 
particles. 

Close to equilibrium, the difference in the rates of scattering into and out of the 
condensate are related by the factor [1 — e'^^'^""'')], where {nc—n) denotes the difference 
between the condensate and the thermal chemical potentials, with determined from 
a time-independent GPE. This gives rise to the following rate equation for the mean 
number of atoms, A^o in thc^ condensate [214, 236] 

^ = 2W+ [(l - e('^=-'')/'=«^) TVo + l] . (210) 

The quantity denotes the scattering rate into the condensate, which, upon 

ignoring all spatial dependence, can be approximated by [214, 130, 234, 236, 237, 240] 

W+ = j dka j dka j dk45(k2 + kg - k4) 

X 6{e2 + £3 - £4 - ^ic)h^2{h + 1) « few x ^""^"^f^^' , (211) 

as briefly explained below. 

Although such a rate can be calculated more precisely, initial studies [213] were 
restricted to somewhat crude approximations: in first instance, apart from ignoring all 
spatial dependence, the non-condensate band distribution /j was approximated by a 
classical Maxwell-Boltzmann distribution e~^^^^~^\ and the integrals were calculated 
over the entire energy range (instead of only within the non-condensate band) . Making 
the additional approximations that fi <C 1, such that fi + lK.1 led to a constant 
rate flxed only by the temperature and the interaction strength, and explicitly given 
by Am{akBT)'^ /irh? , in which the (weak) dependence on the number of atoms has 
been suppressed [214, 236, 243]. Despite the drastic approximations made, this study 
nonetheless revealed good qualitative agreement with experiments. This model was 
subsequently improved to also account for low-lying modes within the condensate 
band, with the full Bose-Einstein distribution used and the range of integration 
restricted to outside the condensate band [237]. An analytical expression obtained 
in the limit when all spatial dependence was ignored yielded a similar qualitative 
description, but with an additional prefactor of few « 3 in Eq. (211). 

In the ergo die approximation, a number of distinct levels is grouped in groups 
of mean energy e/c with Uk = n(efe) = gkf{ek), where denote the energies of the 
dressed levels in the condensate band. In this limit, the occupation numbers evolve 
according to [235, 237] 

^ = [(l - e^^^-^y"-^) + gm] , (212) 

where the rate is now re-expressed as an energy integral over the energies, e^. 

The distribution functions appearing in the integral now become time-dependent, and 
their evolution can be calculated by the ergodic quantum Boltzmann equation [122] 

ijk 
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X 5{ei + e„ - ej - e^) , (213) 

where w is the mean harmonic oscillator frequency of the trap. 

The above equations only deal explicitly with the ^growth! terms, describing 
collisions which lead to a transfer of atoms from one band to the other. Actually, 
there is another process that should be considered here, corresponding to ^scattering' 
terms whereby an atom in the condensate band interacts with another atom in the 
non-condensate band, with population redistribution taking place within each band, 
but without any population transfer between bands [234, 237]. These treatments were 
not explicitly dealt with in previous kinetic formulations. While we do not give their 
respective expressions here (see Appendix C), we note that their presence may lead 
to an acceleration of the initial (spontaneous) growth phase i.e. before growth due to 
bosonic stimulation dominates the dynamics [120]. 

This approach led to the first quantitative predictions of condensate growth [236] , 
and to good agreement with various experiments [237, 239, 244, 238, 245] (see Sec. 
7.2). 

The 'High-temperature Master Equation: ' The exact master equation of Gardiner- 
Zoller for the condensate band (Appendix C.l) may be exploited further, by mapping 
it onto a probabilistic Fokker-Planck equation [215, 216, 217, 241] which is similar, 
but not identical, to that of Stoof. This equation is valid in the regime when 
the eigenfrequencies of the condensate band operator are small compared to the 
temperature, i.e. huo <^ ksT, which roughly coincides with the criterion for large 
occupation per mode. Mapping onto a stochastic differential equation yields an 
equation similar to that of Eq. (200), but with additional contributions due to the 
scattering terms mentioned earlier [215, 216, 217]. By modelling these non-local 
scattering terms approximately [216] (although more accurate expressions can also 
be given [216, 217]), one arrives at a local stochastic GPE for the condensate band 
which takes the form 

da{r, t)= - ^Lca{r, t)dt + Vc [Kair, t)dt + dWair, t)] 

+ Vc[KM{T,t)dt+ia{T)dWM{T,t)] , (214) 

where Vc is a suitable projector into the condensate band, and Lc is essentially the 
unperturbed energy in the condensate band (modified by the mean field of the non- 
condensate). The noise sources appearing here are independent, with dWciY, t) being 
complex and rfVFG(r, t) real. The precise forms of the noise correlations and expressions 
for the redistribution (iCG(r,t)) and scattering (i^M(r,i)) contributions are given 
in Appendix C.2. The formulation of this stochastic equation essentially merges 
the ideas of quantum kinetic theory with those of the Projected GPE formalism, 
whose finite temperature generalisations [119] also include (qualitatively) the coupling 
to a heat bath [216]. The Stochastic GPE has been applied to the dynamics of 
hydrogen condensates [215], to the spontaneous formation of vortices and vortex arrays 
[215, 217, 221] and to equilibrium properties of finite temperature Bose gases [241]. 

At this stage, we wish to point out that, in the limit where the scattering 
contributions are ignored, this equation is identical to that of Stoof (Eq. 200), as 
shown explicitly in Appendix C.2. Note however that additional subtle differences 
do exist between the approaches of Stoof and Gardiner et al. (e.g. determination of 
energy 'cut-off' between the two bands, use of projectors), and the reader is referred 
to the discussion given by Gardiner and Davis [216]. The equivalence between the 
corresponding stochastic GPEs enables us to briefly present a simple damped GPE 
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which arises as a limiting case of both above formalisms [215, 247], and which has been 
used by numerous authors (including one of us - NPP) to study finite temperature 
properties of Bose gases [240, 248, 249, 250, 251, 252]. 

We note here in passing that a related Langevin equation for a single-mode 
condensate was formulated by Graham [246] and used to study fluctuations around 
the equilibrium state of the condensate after it had been formed. 

6.2-4-. The Damped Gross- Pitaevskii Equation: If one heuristically ignores the 
'noise term' r] in Eq. (200), the condensate evolution can be cast (in the classical 
approximation) in the form 

zfi^<i>(r, t) = [1 - *7(r, t, T)] (ho + .g|<I>(r, t)\^ - ^?j <I>(r, t), (215) 

where 7(r, t, T) denotes a dynamical temperature- and position-dependent damping 
rate, given by 7 = i(/3/4)?lS^(r,t) [240, 247]. The addition of a phenomenological 
damping coefficient 7 onto the GPE was originally proposed using general arguments 
by Pitaevskii [248], and first implemented to trapped Bose gases by Choi et al. [249] 
who used a constant, position-independent rate 7 to discuss damping of excitations. 
A similar, but not identical, phcnomenologically-dampcd equation with the factor 
(1 — 27) appearing on the left hand side of the equation has been used in diverse 
studies, including vortex lattice growth [250, 251], and dark soliton decay [252]. 
Although the constant value of 7 chosen in such approaches was such that it agreed 
qualitatively, and to order of magnitude, with experiments, the values used had no 
microscopic justification. The arguments of the preceeding section show how such a 
'phcnomcnologically-dampcd GPE' can be justified from a microscopic perspective, 
while simultaneously providing an explicit expression for this damping coefficient. 

Close to equilibrium, the full self-energy expression of Eq. (199) can be simplified 
to an approximate version which is simpler to deal with. To achieve this, we use the 
identity N{£) + 1 = e'^(^'-''V(e''^^'"''' - 1) = -N{-£), where N{e) = [e^^^^"^) - 1]'^ 
is the Bose-Einstein distribution function. By noting that expressions of the form 
[iV2(A*3-|-l)(iV4 + l)±(7V2 + l)A'^37V4] appear always within integrals containing energy 
and momentum conservation factors (see, e.g. Eqs. (199), (204)), we can immediately 
relate these factors via NiiN^ + 1){N4, + 1)± {N2 + l)N3Ni = [l ± e'^^^-''^ {N2 + 
1)7V37V4 , a result which was used independently by Zaremba, Nikuni and Griffin 
[109], Stoof [220] and Gardiner, Zoller and coworkers [130, 214]. For the purposes 
of our present discussion, we note that, upon considering the expression of Eq. (199) 
for the self-energy sufficiently close to equilibrium, we can use the above relation for 
the sum of the 'in' and 'out' rates, and simultaneously approximate the exponential 
g-/3(e-,i) Rj 1 _ /3(£ _ ^) Rj 1, to obtain 

^ (^'^)"-^7r^ J (2^i (2^y (2^(2.;.) 

X ^ (P2 - P3 - P4) S (Ec £2 - £3 - £4) {N2 + 1)N3N4 . (216) 

For the purposes of the damped GPE it suffices to good approximation to use the 

damping rate 7 = ?'(/'3/4)fiI]^(r, t) w few x 4m,{akBT)'^ /nh'^ deduced by Gardiner, 
Davis and co-workers, where the prefactor has a typical value of around 3 [234, 237]. 

Before concluding this section on alternative beyond mean field approaches, we 
briefly mention two additional methodologies which arc currently receiving increasing 
attention in the ultracold atom physics community, particularly due to their ability 
to handle both weakly- and strongly-interacting regimes. 
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6.3. Positive P- Representation: 

Sec. 6.2 discussed the formulation of a stochastic differential equation in the Wigner 

representation. Alternative choices are also possible, as well-known in the quantum 
optics context [154, 155, 156], and one could for example choose to map the master 
equation into a related stochastic differential equation using the Glauber P function, 
instead of the Wigner function [215]. Although suc;li a function c;annot always 
be interpreted in a probabilistic sense, a variant of this technique, in which the 
number of independent variables is doubled, removes this restriction. Such a positive 
P-representation [253, 254, 255] has the appealing feature that it generates exact 
differential equations, with the absence of explicit vacuum noise contributions in 
the final equations implying that one does not require a mean occupation per mode 
much larger than one, which is a limitation of Wigner- function-based approaches. 
This technique has the advantage that it is exact, and can also handle the strongly- 
correlated regime. The first discussion of this method in the context of BEG was done 
in [195]. Although the validity of its predictions is limited to relatively short times 
(with the sampling error in the simulations growing exponentially after that), there has 
been recent success in studying numerous topics of experimental relevance [256, 257], 
including condensate collisions [258] , with recent developments [259] promising a wider 
applicability in the future. 

6.4. 2-Path Irreducible Effective Action: 

Far from equilibrium dynamics in whic;h initial non-markovian effects become crucial 
can also be described in the context of the 2-Path-Irreducible (2PI) closed time path 
effective action, a non-perturbative technique which has led to significant progress 
in the description of strongly interacting relativistic systems (see, e.g. [260]). This 
approach is also based on a formalism in terms of a Schwinger-Keldysh effective action 
and preserves (at any truncation) important conservation laws such as total particle 
number and energy. Specifically, Rey et al. [261, 262] and Gasenzer et al. [263, 264] 
applied a systematic expansion of the effective action in powers of the inverse numbers 
of field components, which amounts to an expansion of the theory about a strong quasi- 
classical field. This approach, which is also valid for strongly-correlated systems, has 
already been shown to reduce [261] , in the appropriate limits, to the mean field regime, 
the Bogoliubov (one-loop) approximation, and the time-dependent HFB formalism. 

6.5. Brief Summary 

Further approaches have been formulated for the description of ultracold Bose gases at 
finite temperatures which behave either predominantly classically, or are additionally 
affected by quantum effects: 

For the low-lying modes of the system which are typically highly- occupied, one 
can assume that thermal fluctuations largely overwhelm quantum fluctuations, and 
so the system can be described by the (classical) Gross-Pitaevskii Equation. As the 
evolution of this equation is largely independent of the initial conditions ( except for 
very short times), one typically starts from easy-to-implement initial non- equilibrium 
conditions. Fluctuations are evident in single numerical runs, whereas temporal 
averaging (assuming ergodicity) produces smooth density profiles and correlation 
functions. The most detailed implementation also includes a projector to ensure that 
only modes within the 'classical' region under study contribute to the system evolution. 
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In the opposite limit where quantum fluctuations are important, one can still 
pTopagate the system by the same classical equation, while approximately accounting for 
quantum effects by the inclusion of random fluctuations of, on average, half a particle 
per mode in the initial conditions, in the so-called Truncated Wigner approximation. 
At higher tern,peratures, this approach leads to spurious damping, and its validity is 
therefore restricted to relatively short times. 

In general, one can envisage separating the modes of the trapped gas into low- 
lying modes that should be treated accurately, and higher-lying (thermal) modes that 
can be treated semi-classically. The low-lying modes are found to obey a so-called 
Stochastic Gross-Pitaevskii equation; this is a dissipative Gross-Pitaevskii equation 
where dissipation arises from the coupling of these m,odes to the reservoir of high-lying 
modes; unlike the semiclassical treatments discussed earlier, this equation additionally 
includes a dynamical noise contribution with gaussian correlations to account for fluc- 
tuations. Simultaneously, the high-lying modes are described by a Quantum Boltzmann 
Equation. On the formulation side, this latter approach constitutes the appropriate 
generalization of the self- consistent Gross-Pitaevskii-Boltzmann, or 'ZNG' formalism 
- note however that the stochastic Gross-Pitaevskii equation has so far only been nu- 
merically implemented when coupled to a static heat bath, with the implementation of 
the required coupling to the Quantum Boltzmann Equation currently being pursued. 
In the context of the stochastic Gross-Pitaevskii equation, single-run results indicate 
typical fluctuations present in experiments, while suitable averaging over a number of 
runs is required to generate smooth density profiles and correlation functions. 

This concludes our formal presentation of both mean field and more advanced 
formalisms used to describe ultracold Bose gases at finite temperatures, and the 
next section puts these theories into context by briefly comparing their respective 
predictions for experimentally relevant cases. 

7. Applications and Comparison of Above Approaches 

The theoretical study of ultracold atoms has been enhanced by the possibility of 

immediate comparison to experiments, and the theories presented above have been 
applied to diverse experimental conditions, as mentioned throughout this Tutorial. 
Initial applications of these theories focused on static properties, such as condensate 
fractions and density profiles, which arc described fairly well already at the Hartree- 
Fock level [12], with relatively small (but measurable) differences provided one is far 
from the critical region in 3D, or the corresponding regimes of large phase fluctuations 
in ID and 2D. As there arc unfortunately no 'benchmark' tests available to date in 
the ultracold gas community (despite recent efforts by the authors), we restrict our 
comparison below to two key experimental issues which both stimulated and assisted 
in the development of the above approaches; we do not give any new results here, but 
merely reproduce previously published results. The subsequent discussion is therefore 
not intended as a detailed overview of those research topics, but rather as a brief 
introduction to highlight some key issues mentioned in our preceeding discussion; 
accordingly, the list of references given below is far from complete. 
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7.1. Finite Temperature Excitation Spectrum 

Studies of the response of a system to external perturbations has been the subject 
of detailed investigation in 'traditional' condensed matter systems over the past 
decades [265]. The excitations can be split into 'coUisional' (or hydrodynamic) and 
'coUisionsless' regimes, depending on the density of the system, or equivalently on 
the coUisional mean free path: Most of the (early) experiments with ultrac;old gases 
probed the latter regime, in which collisions between thermal atoms must be explicitly 
considered, and the dominant effects arise from self-consistent mean fields. These can 
be further classified depending on the relative size of the excitation wavelength Aox 
compared to the healing length ^ oc (/um)~^ of the system: for Aex ^ (i.e. large 
momenta) one obtains 'single-particle' excitations, whereas in the opposite regime 
of small momenta (Aqx ^ (,) one obtains phonon-likc, or collective excitations (sec 
Sec. 2.2 and Eq. (43)). In confined systems, the size of the trap sets an additional 
lengthscale to the problem [61]. When Aex approaches the condensate size, the 
excitation spectrum becomes discretizcd, i.e. the low energy collective modes of the 
system are standing sound waves at specific frequencies, whereas excitations with 
energies larger than the typical trap frequency behave semi-classically [59]. 

The study of discrete collective excitations provides a very stringent precision 
test for the validity of theoretical approaches. In the early years of atomic BEC 
experiments, most of the measurements (density profiles, expansion dynamics) , could 
be reasonably explained by simple mean field theory [12, 56], with the collective 
oscillation frequencies at low temperatures measured at JILA [266] and MIT [267] well- 
explained by the zero-temperature Bogoliubov equations (Eq. (40)) [56, 60]. However, 
subsequent measurements performed in the presence of a thermal cloud [268, 269] 
posed significant challenges to theorists. In particular, an experiment at JILA [268] 
which measured the temperature dependence of the oscillation frequencies of a slightly- 
elongated condensate was not in agreement with existing theoretical models, and the 
theoretical effort to understand the physics of this experiment was an important drive 
for the establishment of dynamical finite temperature theories. 

The JILA experiment, in which a time-dependent sinusoidal perturbation 
distorted the trap transversally, measured the quadrupole oscillation modes of the 
condensate shown schematic in Fig. 13; these are characterized by the projection m of 
the angular momentum onto the weakly-confining trap axis as follows: (i) in the m = 
cylindrically-symmetric mode, axial and radial directions oscillate out of phase, and 
(ii) in the m = 2 mode the cylindrical symmetry is broken, with the condensate radial 
oscillations out of phase with each other. The oscillation frequencies of both modes 
were observed to decrease with increasing temperature, with the m = mode however 
additionally displaying an 'anomalous' upward shift at temperatures T > 0.7Tc, with 
the observed frequencies differing by about 10 — 20% from the corresponding zero- 
temperature results [268]. The experimental data of the m = mode, along with 
predictions based on various theories applied to this particular problem are shown in 
Fig. 14, and arc explained in detail below. 

The first finite temperature attempts to model this experiment were based on 
a static thermal cloud. Initially, the static HFB theory was applied in the so-called 
'Popov' limit of no anomalous correlations (Sec. 3.3.1); this was found to provide an 
excellent description for temperatures up to 0.6Tc (Fig. 14, '-|-' symbols); remarkably, 
the HFB-Popov predictions were found to be only slightly different from corresponding 
T = results when one additionally accounted for the different number of condensate 
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SCHEMATIC OF EXCITATION FREQUENCIES MEASURED AT JILA 

Quadrupolar Oscillation Quadrupolar Oscillation 

m = I m 1= 2 




Axial & radial lengths of condensate Radial lengtlis of condensate 

oscillate out of phase oscillate out of phase with each other 



Figure 13. (colour online) Schematic of two different quadrupolar modes of 
excitations observed in early BEC experiments [266, 267, 268, 269], with emphasis 
here on the modes observed in the slightly elongated condensates at JILA. Each 
excitation is identified by the angular momentum of the excitation about the z- 
axis, denoted by m. The m = mode corresponds to a cylindrically-symmetric 
quadrupolar excitation, whereas the cylindrical symmetry is broken for |m| = 2, 
with axial and radial motion taking place out of phase with each other. While 
the m = mode was well understood, experimental attempts to excite the m = 2 
mode at finite temperatures generated some 'anomalous' features discussed in the 
text. (Image modified from similar picture presented in [61]). 



atoms at the corresponding temperatures [270, 271]. Tlie need to understand tlie 
system behaviour at higher temperatures stimulated the development of improved 
theoretical models: motivated by analytical work performed also by one of us (NPP) 
[85], Hutchinson, Dodd and Burnett carried out simulations which explicitly included 
the anomalous average and approximate many-body effects [90] via the so-called 
generalized HFB models (Sec. 3.4). From the two models proposed, the one based on 
replacing the two-body effective interaction strength g by g{v) throughout Eqs. (82)- 
(85) (i.e. both for condensate-condensate and condensate-thermal collisions) led to a 
correct prediction of the downward shift of the m = 2 mode; nonetheless, such a theory 
incorrectly predicted a similar shift for the m = mode (Fig. 14, '*' symbols). Similar 
results were obtained by Minguzzi and Tosi based on a time-dependent linearized 
Hartree-Fock model [272], by Shi and Zheng in the context of a finite temperature 
variational method [273] and by Reidl et al. [274] by means of the dielectric formalism. 
As already mentioned, these models had the common feature of a static thermal cloud, 
which was quickly realized to be too restrictive for the particular experiments. In 
order to overcome this limitation, Giorgini performed a linear response treatment of 
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Reduced Temperature T / T 



Figure 14. (colour online) Comparison of the predictions for the temperature 
dependence of the excitation frequency of the m = mode measured at JILA 
[268] (shown by black squares) based on different theoretical models. Plotted 
figure is restricted to the regime T > O.STc, with the different theories at 
lower temperatures generally converging to the experimentally observed (and 
theoretically predicted) T = value of 1.86. Blue (filled): Number-conserving 
formalism with (circles) or without (diamonds) direct excitation of the thermal 
cloud from the probe [167]. Green: Static thermal cloud theories with ('*', 
generalized HFB with gt = g(r) [90]) or without ('+', HFB-Popov [271]) inclusion 
of the anomalous average. Red (filled): Predictions of ZNG approach for different 
excitation probe frequencies aimed at exciting primarily the condensate (circles, 
LD = 1.75ij;z) or the thermal cloud (inverted triangles, uj = 2u)z) [136]. Brown: 
In-phase (left triangles, top) and out of phase (right triangles, bottom) mode of 
excitation between condensate and thermal cloud [275, 276]. (Note that all data 
have been extracted manually from corresponding published works; each data 
point has an independent small error which should be on the order of the size of 
the symbol). 



the coupled dynamics (Sec. 4.3.1) but nonetlieless found a similar downward shift for 
both modes [100]. 

A qualitative explanation of the observed 'anomalous' temperature dependence 
of the m = mode was given by Bijlsma, Al Khawaja and Stoof in the context of a 
finite temperature GPE in the Hartree-Fock approximation coupled to a coUisionless 
Quantum Boltzmann equation [275] , a study subsequently generalized to the coUisional 
regime [276]; although these works were motivated by the full theory of Stoof 
(Sec. 6.2.2), the analysis was performed here in the simplified semi-classical mean 
field limit which essentially amounts to the 'ZNG' theory presented in Sec. 4.4.3. 
However, instead of solving these equations self-consistently, their analysis was 
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based on a variational approach. Their conclusion was that the reported m = 

excitation frequency actually corresponded to simultaneous measurements of two 
distinct excitations (brown left/right triangles in Fig. 14) which become coupled 
above a certain temperature; a similar interpretation was also given by Olshanii [277] . 
Bijlsma, Al Khawaja and Stoof suggested that these modes corresponded to in-phase 
and out-of-phase oscillations of the condensate and the thermal cloud. This picture was 
later confirmed by one of us (B J) by detailed numerical simulations based on the 'ZNG' 
theory (Sec. 4.4.3) which accounts for the full dynamical coupling between condensate 
and thermal cloud; this work identified the two distinct modes in which the condensate 
is excited as corresponding to a damped condensate oscillation in contact with a quasi- 
static heat bath, and an oscillation coupled to the motion of the thermal cloud [136]. 
This work highlighted the extreme sensitivity of the observed system response to the 
precise details of the imposed perturbation (filled red inverted triangles /circles in Fig. 
14), such as the driving frequency which was not very accurately determined in the 
experiments. The temperature dependence of the damping rates based this theory 
was also found to closely match the experimental observations. The magnitude of 
the observed damping rates was also interpreted by means of a perturbative approach 
[101, 102] and field theoretic techniques [278]. 

The full proof of the simultaneous excitation of the two modes in the experiment, 
was later presented by Morgan, Rusch, Hutchinson and Burnett [166, 167, 168] 
by an extension of their earlier work [78, 279] which accounts for all second order 
processes within a number-conserving formalism (Sec. 5.3). They explicitly showed 
(blue circles/diamonds joined with dotted lines in Fig. 14)) that the 'anomalous' 
behaviour of the m = mode is a result of a temperature-dependent interplay in 
the method of excitation of the condensate: at low temperatures, this is dominated by 
direct excitation of the condensate from the probe, whereas at higher temperatures 
the condensate is excited predominantly by its coupling to the thermal cloud, which 
is itself excited from the probe. This detailed study led further to the identification 
of novel dynamical resonances which could be experimentally observed [169]. 

In addition to the above comparison to experiments, another crucial test of the 
consistency of a kinetic theory, is provided by the so-called Kohn mode. Kohn's 
theorem for electrons states that the cyclotron frequency is not affected by interactions 
in a static magnetic field [280, 281]. By analogy, the dipole mode of the condensate, 
which corresponds to the centre of mass oscillation, should occur at the harmonic 
oscillator frequency, without being affected by 2-body interactions, since in a harmonic 
trap the centre of mass is explicitly decoupled from the internal degrees of freedom. It 
turns out that the Kohn mode is only accurately reproduced by theories which treat 
the thermal cloud dynamics on an equal footing as that of the condensate. These 
include the second order excitation theory of Morgan and Burnett [166, 167, 168], 
the 'self-consistent Gross-Pitaevskii-Boltzmann' or 'ZNG' approach [109] and the fully 
dynamical stochastic approaches of Stoof [128] and Gardiner-Zoller-Davis-Ballagh and 
co-workers [130, 216]. 

Our preceeding discussion was restricted to a specific mode of excitation, as a large 
number of diverse theories was applied to its study, thus enabling a direct comparison 
between them. Numerous other excitation frequencies were studied in detail using 
diverse theories, and we briefly remark here that, where applied, the self-consistent 
Gross-Pitaevskii-Boltzmann, or semi-classical 'ZNG' theory of Sec. 4.4.3 has yielded 
excellent agreement with such experiments [132, 134]. 
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7.2. Condensate Growth 

The study of quantum phase transitions has been an active and challenging research 
topic in diverse areas of physics. In the context of ultracold Bos(^ gases one is interested 
in understanding the process of condensate formation and the associated evolution of 
coherence, leading to the establishment of off-diagonal long-range order. Apart from 
being an interesting topic from a fundamental point of view, suc;h studies pose the 
most stringent validity test of the various theoretical models of ultracold Bose gases. 
It is thus appropriate to discuss here to what extent the presented models facilitate 
an understanding of the strongly non-equilibrium features obsc^rved experimentally. 
Early theoretical work on this issue (initiated even before the pioneering experiments 
in 1995) focused mainly on (i) a qualitative identification and understanding of the 
distinct stages of the process of condensate formation and equilibration starting from 
a thermal gas [127, 172, 173, 174, 175, 212, 282, 283, 284, 285, 286, 287, 288], and, 
more specifically (ii) on the evaporative cooling process [121, 122, 289, 290, 291, 292] 
which plays a dominant role in the experimental route towards condensation. Clearly 
the topic of condesate formation is challenging from a theoretical perspective, and only 
the most advanced approaches presented in this Tutorial may be realistic candidates 
for succeeding in such a task. 

Pioneering work on this issue was carried out in 1998 at MIT, where in situ 
studies of condensate formation were performed [244]. This experiment introduced the 
technique of 'shock cooling', whereby a trapped thermal cloud is suddenly quenched 
below the transition point, and the evolution of the resulting 'supersaturated cloud' 
into a condensate is studied. Theoretically, this initial condition can be modelled by 
truncating the Bose distribution at some temperature T, corresponding to an energy 
ksT above which all atoms are assumed to be efficiently removed by evaporative 
cooling. The condensate was observed to grow slowly initially (spontaneous growth), 
before the bosonic enhancement set in, leading to exponential growth. This process 
was shown to be distinct from a simple relaxation mechanism, with the condensate 
eventually equilibrating when its chemical potential approached that of the thermal 
cloud [244]. 

Two theoretical approaches have been applied to model this experiment, as 
discussed below, although at this point we should also highlight the related work 
of Barci, Fraga, Gleiser and Ramos on the growth and equilibration of a homegeneous 
Bose condensate [230, 231, 232]. The first numerical studies of condensate growth 
came from Gardiner, Zoller, Ballagh and Davis [236], before any specific condensate 
formation data were available. Their analysis was performed in the context of the 
'simple' quantum-optical rate equation (Eq. (210)), and appeared to be in qualitative 
agreement with existing experiments [1, 2]. Following the pioneering MIT controlled 
growth experiment [244], their analysis was improved to additionally include the 
dynamics of the occupations of low-lying trap levels via Eqs. (212)-(213); this provided 
the first quantitative results which yielded good overall agreement with the experiment 
[237]; furthermore, treatment of the scattering terms was found to speed up the 
initial period of growth. Their model was further improved by the inclusion of 
temporal depletion of the non-condensate [235] , which effectively modified the simple 
rate equation by the replacement of the static thermal cloud chemical potential 
by a dynamical effective chemical potential determined self-consistently, leading to 
improved agreement [235]. 

Parallel work was undertaken by Bijlsma, Zaremba and Stoof in the context of 
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Figure 15. Compajrison of a set of experimental data from the MIT group 
(filled circles) [244] and theory. Shown are apparent fits to the data generated by 
Bijlsma, Zaremba and Stoof in the context of the self-consistent Gross-Pitaevskii- 
Boltzmann or 'ZNG' theory (dashed line), and by Davis, Gardiner and Ballagh 
based on the condensate master equation (solid grey); both these simulations - 
which are practically identical except for small differences in the initial growth 
characteristics due to the different initial condensate numbers assumed in each 
implementation - are performed for an initial thermal cloud of 40 X lO'' •^'^Na 
atoms at an initial temperature of Ti = 765nK, and a relatively severe truncation 
of the energy distribution. Equilibrium predictions resulting from these initial 
parameters appear to be in disagreement with the observed static (equilibrated) 
MIT data. Choosing instead the parameters so as to match the static experimental 
data and performing simulations with the same condensate master equation 
actually leads to very different growth dynamics, as shown by the solid black 
line. (Reprinted figure with permission from M.J. Davis, C.W. Gardiner and R.J. 
Ballagh, Phys. Rev. A 62, 063608 (2000). Copyright (2000) by the American 
Physical Society.) 



the 'self-consistent Gross-Pitaevskii-Boltzmann' or 'ZNG' theory (which constitutes a 
special case of the theory put forward by Stoof in Sec. 6.2.2); these authors arrived at 
similar conclusions [131], with their predictions found to be compatible with the results 
of Davis, Gardiner and Ballagh for approximately the same initial parameters [235]. 
Before comparing the predictions of these two theories to the experiments, it is worth 
remarking that although the semi-classical mean-field model cannot describe the onset 
of condensation, it can nonetheless accurately describe the subsequent dynamics; thus 
the semi-classical theory can still provide good agreement with experiments provided 
the condensate is initially artificially seeded with a small number of atoms (typically 
determined by the number of atoms in the lowest harmonic oscillator state at the 
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temperature of the truncated Bose distribution [131]). 

The predictions of these theories are compared against each other and a particular 
set of experimental data from the MIT group (corresponding to Fig. 4 in [244]) 
in Fig. 15. Although the agreement between these two numerical implementations 
(denoted respectively by solid grey and black dashed lines) was found to be very good 
[235], and even appears to reproduce the experimental data (dots) rather accurately, 
we note that the displayed simulations are not actually for parameters consistent 
with the experimental data quoted by the MIT group (or inferred from their data). 
Here, it should be noted that, in theoretically analyzing the experiment, it was not 
always clear what input parameters the model should use in order to correspond to 
the experimental conditions, as not all such parameters had been measured in (or 
extracted correctly from [234]) the experiment. This creates certain ambiguity in 
performing detailed comparisons. While most of the experimentally obtained curves 
seemed to be generated by ab initio application of the theory (using 'admissible' 
experimental parameters), there was at least one 'puzzling' case (shown in the 
figure) where the theoretical model could not reproduce the observed behaviour for 
parameters consistent with the experimental values [235]. In particular, a simulation 
performed by Davis, Gardiner and Ballagh for parameters extracted from the MIT 
static experimental data (after equilibration) actually generated the solid black curve 
(instead of the solid grey one) - in clear disagreement to the measured data. Despite 
the good qualitative and even in some cases quantitative agreement, the fact that 
not all observed growth curves could be reproduced ab initio from the experimental 
parameters, required further study, such as removing the condition of ergodicity (as 
undertaken by one of us - BJ - in collaboration with Zaremba) [137], incorporating 
the quasiparticle nature of the spectrum of low-lying modes, or allowing for a non- 
adiabatic condensate growth exhibiting shape oscillations. 

The next systematic experimental study of condensate growth was performed in 
2002 by Kohl et al. for a ^''Rb condensate [238]. This study was different from the 
MIT one, in that it was performed under continuous evaporative cooling into the 
quantum degenerate regime, thus enabling (for suitable parameters) to observe the 
formation process in 'slow motion'. In this work, all relevant parameters required 
for an ab initio modelling of the experiment were measured. In the case of strong 
cooling (as in the MIT experiment), this led to remarkable agreement with the theory 
of Davis, Gardiner and Ballagh with no free parameters (once the modification of 
the trap potential due to gravity was taken into account [239]). However, Kohl et 
al. also studied the opposite regime of slow cooling., in which the system spends a 
longer time around the transition region, and found a 'two-stage' growth [238]: the 
evolution featured a slow initial growth which could not be explained by the above 
model, with this regime followed at later times by the usual exponential growth to 
equilibrium. As these features occured very close to T^, it was suggested that the 
slowing down of the initial growth rate could arise as a result of the appearance of a 
quasi-condensate, which is thought to preempt condensate formation even in three- 
dimensional systems [172, 173, 175, 286]. Subsequent related experiments performed 
at Amsterdam [293], Orsay [245] and Zurich [294] provided contradictory indications 
regarding the appearance and role of the quasi-condensate in such systems, which is 
believed to be associated with the appearance of condensate shape oscillations. 

Such shape oscillations were indeed observed, and strong phase fluctuations mea- 
sured in an experiment studying the formation of condensation into nonequilibrium 
states at Amsterdam, which highlighted the role of local thermalization [293] . The Or- 
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say experiment measured the evolution of both population and coherence in a highly 

elongated trap. While they found evidence of shape oscillations, there was no indica- 
tion of a two-stage growth despite being deep within the quasi-condensation regime 
[245]. In this limit, theoretical models which do not fully account for phase fluctuations 
cannot a priori be assumed to be valid. Nonetheless, (perhaps somewhat unexpect- 
edly) the model of Eqs. (210)-(213) led to good agreement with experiments, up to an 
arbitrary delay time in the onset of condensate growth [245] ; while such a delay time 
was introduced in an ad hoc manner in order to match the observed growth curves, 
it is still unclear whether this delay is related to phase fluctuations. Moreover, recent 
work in the real-time observation of the formation of long-range order performed in 
Esslinger's group found no evidence of a quasi-condensate stage [294]. Thus, there ap- 
pears to be some controversy regarding the appearance and effect of shape oscillations 
and phase fluctuations in the evolution of a quenched strongly-non-equilibrium Bose 
gas into an equilibrium condensate. This and related issues require more advanced 
studies and may be resolved by simulations based on the stochastic techniques of Sees. 
6.2.2-6.2.3. 

Simulations of the stochastic GPE of Stoof (Eq. (200)) have in fact been performed 
in collaboration with one of us (NPP) in the context of condensate growth in a dimple 
microtrap [223] (and in the related study of atom laser dynamics [225]); in such 
experiments, the onset of condensation is induced by an entirely different mechanism, 
based on adiabatic local phase-space compression [295] . The first such experiment was 
performed at MIT [226]: The harmonic confinement of a thermal cloud at equilibrium 
was perturbed by the addition of a narrow gaussian dimple trap at the centre of 
the cloud, yielding a maximum increase in phase space density by a factor of 50. 
Condensate formation in the limit of slow trap addition was studied by means of non- 
destructive imaging, revealing condensate fractions of up to 20% for different dimple 
depths. Importantly, application of a sinusoidal modulation of the trap depth led to 
a controlled and reversible crossing through the phase transition, with the condensate 
formation found to lag the trap modulation by about 70 ms [226]. 

The very first application of the stochastic GPE of Eq. (200) to ultracold gases 
was in fact undertaken by Bijlsma and Stoof [220] in a purely one-dimensional context 
in an attempt to analyze this experiment; this yielded good qualitative agreement with 
the experimental results with respect to the magnitude of the previously mentioned 
lagging time. However, a detailed quantitative comparison could not be performed, 
as the experiment was undertaken in a fully-three-dimensional regime, whereas the 
simulations were limited to the one-dimensional regime. It is important to remark 
here that attempts to model the dynamics by a semi-classical theory consisting of 
a finite temperature GPE dynamically coupled to a reservoir of thermal atoms, i.e. 
the stochastic GPE of Eq. (200) without any noise, failed to describe the experiment 
accurately; in particular, such simulations predicted cycles of successively decreasing 
numbers of condensate atoms, with the rate of decrease of the condensate atom number 
between such cycles depending on the initial conditions. The above scenario highlights 
one characteristic example where stochastic theories are required in order to accurately 
describe the observed behaviour. 

In the same experimental publication, Stamper-Kurn et al. [226] suggested 
studying non-equilibrium (quasi) condensate growth induced by the sudden addition 
of a dimple trap. Preliminary experimental work along this line was undertaken 
in the group of Jorg Schmiedmayer by the addition of a deep dimple microtrap on 
the weakly-confining axis of an atomic gas trapped close to an atom chip [296]. 
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Theoretical modelling of this scenario was performed by one of us (NPP) in the 
context of the stochastic GPE (Eq. (200)), focusing on the case of a weakly-interacting 
quasi-one-dimensional Bose gas [223]. This led to the identification of a range of 
competing dynamical phenomena (shock wave propagation, direct quasi-condensate 
growth) which could be detected in future experiments. 

Other characteristic dynamical examples where beyond mean field theories are 
required include (but are not restricted to) (i) the spontaneous process of vortex 
formation (see Fig. 12) studied recently in more detail in [221], and the related issue 
of the dynamics of the Berezinskii-Kosterlitz-Thouless (BKT) transition [80] in two 
dimensions to a state which does not possess long-range order but is instead associated 
with the unbinding of vortex pairs [194]. 

While it is believed that the stochastic approaches, along with the appropriate 
time-dependent treatment of the thermal cloud should successfully address these 
issues, this presumption will need to stand the test of time, particularly as 
experimentalists are continuously exploring novel interesting phenomena which 
generally provide more stringent tests to the existing theories. 

8. Conclusions and Outlook 

This Tutorial has presented a desc;ription of the most common theoretical methods 
available for describing weakly-interacting ultracold Bose gases at finite temperatures, 
and some concluding remarks are needed here. 

Before doing so however, we feel compelled to note that, although there is 
consensus on a limited number of points (e.g. the evident importance of the thermal 
cloud, or the need to account for phase fluctuations in low dimensions), there is 
no universally accepted 'optimal' theory for the description of such systems, with 
researchers coming from different communities typically taking different 'stands' on 
this topic. This makes such a concluding discussion hard and even possibly somewhat 
'controversial'. In our preceeding presentation we have gone to great length to ensure 
objectivity of discussion and avoid introducing any personal bias into this Tutorial, 
and we hope our conclusions reflect the same impartiality. To aid the less experienced 
reader, we have attempted here a rough classification in terms of the Table shown in 
Fig. 16. Clearly such a presentation, while indicative, should not be taken too literally, 
as it is impossible to summarize an entire theoretical development in a single sentence, 
let alone classify each theory by 'ticking relevant boxes'. With these points in mind, we 
proceed with some concluding remarks, all of which relate to dilute weakly-interacting 
gases with na^ <C 1 and a <C XdB, such that only elastic binary collisions are relevant 
for describing the system properties. 

8.1. Classification of Theoretical Approaches: Role of Symmetry-breaking 

The first point to make is that, loosely speaking, one can classify existing theoretical 
formalisms into three different 'classes' of approaches, based on certain common 
conceptual notions shared between them. 

The main such notion is the idea of symmetry-breaking: In both high-energy 
physics and condensed matter physics it is convenient to assume that the operator 
describing a relevant quantum field reduces, below the critical temperature for some 
phase transition, to an appropriate macroscopic complex (wave)function. This implies 
that the macroscopic variable no longer exhibits the symmetry (or all symmetries) 
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of the original system hamiltonian, which (for the systems we are interested in) is 
invariant under a change of phase of the operator; in this case one says that the 
symmetry has been (spontaneously) broken. While views are generally divided as to 
whether this is a 'physical reality' or simply a convenient 'mathematical concept', in 
some systems this makes no difference as there is no definite (fixed) number of particles 
- that is certainly the case in high energy physics, and the same is essentially true for 
conventional (macroscopic) condensed matter systems. In such cases, the concept of 
symmetry breaking has proven extremely useful in explaining the system's behaviour. 
The notion of broken symmetry however becomes somewhat more 'problematic' in the 
case of a small (finite) system of atoms, as is relevant for experiments with ultracold 
trapped gases, which typically consist of no more than a few (tens of) millions of 
atoms; the reason is that if one were to accept the view that at any time the system 
consists of a definite number of atoms, then the mean value of the Bose field operator 
^'(r, t) should be identically zero. 

The simplest - what one might call generic - approaches for ultracold gases rely 
on symmetry-breaking, i.e. on the assumption that the Bose field operator ^'(r,t) can 
be split into a mean-field condensate contribution (/)(r, t) (denoting the condensate 
wavefunction) and an operator describing the (quantum and thermal) fluctuations 
about this mean field. Approaches relying on this separation are known as 'mean- 
field' approaches, whereas treatments which explicitly maintain the operator nature 
of the condensate contribution (i.e. 0(r,t) (j){r,t)) automatically conserve the 
total atom number, and are therefore known as 'number-conserving'. In a large 
coherent three-dimensional system far from the regime of critical fluctuations such 
a distinction becomes an 'academic issue', i.e. one is essentially 'free to choose' the 
'class' which one is most comfortable with, with predictions between the two classes of 
approaches believed to lead to negligible differences (which may not even be observable 
experimentally). The same is however not true in all experimentally-relevant regimes. 

8.2. Mean Field Theories 

Mean field theories are typically based on identifying a set of slowly- varying mean- field 
quantities and formulating a closed system of equations to be solved self-consistently 
either in a static or a dynamical manner. However, such a closed system is not 
exact, and to arrive at such a formulation one has to impose certain decorrelation 
approximations, i.e. assume that averages of products of operators denoting the 
fluctuations about this mean field can be split into products of 'lower-order' averages, 
each containing a smaller number of operators. Depending on the 'crudeness' of the 
imposed approximations and the related choice of the 'mean field variables' of the 
system, one can generate all different mean field theories, as discussed in Sees. 2-4, with 
such treatments essentially formulated in first or second order perturbation theory. 
First order perturbation theory corresponds to interactions of static self-consistent 
mean fields, whereas formulations to second order in the effective interaction strength 
also allow for dynamical (particle) exchange between such generalized mean fields. 
Note that second order perturbation theory works well here when combined with a 
pseudopotential approximation to the upgraded effective interaction (T-matrix), as 
it can essentially be interpreted as an expansion in terms of the difference of the 
actual effective interaction experienced by two atoms in the medium of condensed and 
thermal atoms from the corresponding one (two-body or many-body T-matrix) used 
within the particular theory under consideration. 
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Mean field theories can be typically grouped into different categories, depending 
on the following three 'classification issues', with each of these corresponding to 
different physical content: 

(i) Nature of Excitation Spectrum: Are the system excitations treated as single- 
particle ones (modified by mean fields to the level of approximation), or does 
one need to consider 'particle mixing' into quasiparticles, which exhibit a different 
(phononic) excitation spectrum at low momenta? Inclusion of quasiparticle physics is 
a straightforward, yet lengthy process, and is lacking in most current implementations; 
nonetheless, this does not necessarily place huge restrictions on the validity of these 
theories for dilute weakly-interacting gases, unless they are applied to systems at 
extremely low temperatures In the finite temperature theories of interest, this choice 
respectively makes a distinction into 'Hartree-Fock'- and 'Hartree-Fock-Bogoliubov'- 
type of approaches. 

(ii) Anomalous Average and Many-Body Effects: At finite temperatures, one 
should at least describe the system by a mean field for the condensate and a 
corresponding one for the non-condensate, or thermal cloud. If one additionally 
chooses to include the so-called (pair) anomalous average as a 'fundamental' system 
variable, the theory will then incorporate additional physical content, in the sense 
that many-body effects, which lead to a modification of the effective scattering length 
due to the presence of the medium, are largely accounted for. Note however that 
extreme care is needed when doing so, as not all resulting theories treat atom-atom 
interactions consistently, and may even lead to worse predictions than 'more basic' 
theories for certain quantities of interest. 

(iii) Particle- exchanging Collisions: Most importantly, does the theory only 
include mean-field coupling between the condensate and the thermal cloud 
(formulation to lowest order in the effective interaction), or does it additionally allow 
for the dynamical exchange of particles between these two sub-systems (second order 
formulation)? In the former case of purely mean field coupling, the theory is limited 
essentially to equilibrium properties, e.g. density profiles and condensate fraction: 
there is a diverse range of theories to choose from here, depending on the issues 
raised earlier, with possible choices (starting from the most basic) being Hartree- 
Fock, Hartree-Fock-Bogoliubov(Popov), or generalized Hartree-Fock-Bogoliubov. If 
one is only interested in bulk quantities or qualitative features, then Hartree-Fock 
is an excellent starting point - although a detailed comparison to experiments will 
typically require (at least) use of the other more generalized mean field theories. In the 
other extreme, i.e. when one is interested in quantities which depend critically on the 
dynamics of either (or both) of the sub-systems, one should necessarily include atom 
transfer collisions between the condensate and the thermal atoms, and, in general, 
also redistribution collisions within the thermal cloud. In this case, one should best 
resort to a description in terms of a dissipative Gross-Pitaevskii equation which is self- 
consistently coupled to a 'reservoir' of high-lying modes described by the Quantum 
Boltzmann equation; the latter, is an extension of the usual Boltzmann equation for 
the distribution function of a gas when it is dynamically coupled to a condensate, and 
is somtimes referred to as the 'ZNG' theory. In addition to constituting a suitable 
choice as long as one is not too close to the regime of critical fluctuations, such a 
theoretical description also enables one to 'switch off' at will either, or both, of these 
coUisional contributions to investigate their relative role. Such investigations have 
been undertaken in various studies (also by the authors - mainly BJ), and in general 
the mechanism of transfer of atoms between the condensate and the thermal cloud is 
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crucial in order to accurately predict experiments looking at dynamical effects, such 
as damping of elementary or macroscopic excitations. 

In order to derive all mean field theories mentioned above, one essentially 
separates the full system hamiltonian accounting for elastic binary interactions into 
contributions depending on the number of non-condensate operators (from zero to 
four) appearing within each contribution. Terms up to quadratic order only accurately 
take account of the condensate mean field and small excitations on top of this, and 
are generally limited to T w 0; treating the remaining higher order contributions 
purely within 'suitable' mean field approximations leads to an approximate inclusion 
of thermal effects, but only within the context of static variables, i.e. self-consistent 
mean field coupling. Inclusion of particle-exchange collisions additionally requires the 
perturbative treatment of the difference between the exact higher-order contributions 
to the hamiltonian from their corresponding approximate expressions based on the 
selected mean field approximations. 

However, it is important to note that no mean field theory is able to predict the 
behaviour of the system at the critical point, or to simulate condensate growth from 
a purely thermal initial gas (although a self-consistent Gross-Pitaevskii-Boltzmann or 
'ZNG' model nonetheless works rather accurately as soon as there is a small condensate 
'seed' present in the initial conditions). In this case, one requires consideration 
of beyond-mean-field, or number-conserving approaches, in which the condensate 
contribution is described by an appropriate operator, rather than by a corresponding 
mean field. There are essentially two 'types' of such approaches, as highlighted below: 

8.3. Number-conserving Perturbative Treatments 

The first such approach is very similar in procedural development to the treatment of 
individual contributions to the main system hamiltonian mentioned above, although 
there are important conceptual differences. The main difference here is that the 
excitations are by construction orthogonal to the condensate (not strictly true in 
mean field approaches), and the non-condensate operator is defined in a slightly 
different manner which explicitly conserves the total number of atoms. Apart from 
not requiring the use of symmetry-breaking (i.e. the condensate contribution is still an 
operator), at the expense of reasonable algebraic complexity, such treatments are well- 
suited for small trapped gases where such effects and issues associated with precise 
number-conservation may become important. The power of such approaches has been 
demonstrated in distinguishing between direct and indirect excitation mechanisms for 
condensate oscillations in a particular experiment; however, such a formalism has 
not yet been widely applied and (in light of the significant extra effort involved in 
formulating/simulating the theory) its benefit over other approaches remains unclear; 
this issue will presumably be resolved when such approaches are generalized to the full 
dynamical inclusion of the thermal cloud within the context of a fully self-consistent 
second order perturbation theory which is, however, presently lacking. 

8.4- Alternative Number-conserving Approaches: Classical Field and Stochastic 

This brings us to the third class of approaches, which incorporates classical-field and 
stochastic methods. A key difference of these to all earlier approaches is that such 
treatments include fluctuations (in the form of added numerical noise for the stochastic 
approaches), and the results for various parameters should be extracted via suitable 
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averaging - either temporally, or over a number of different numerical realizations. 
In addition, these theories describe the low-lying modes of the system in a 'unified' 
manner, so the identification of the 'condensate', or 'quasi-condensate' in these theories 
is slightly different to the perturbative mean field, or number-conserving approaches 
discussed previously - an important point to bear in mind when attempting direct 
comparisons between different theoretical formulations. 

There are three main approaches)} we wish to highlight here: 

(i) The Classical Field Method: In this approach, as suggested by its name, the 
entire evolution of the system is described classically. Such a treatment should be valid 
as long as thermal effects become so large, that they essentially 'wipe out' all quantum 
corrections. The classical field describing the system is propagated via the Gross- 
Pitaevskii Equation (GPE), with the possible addition of a projector to ensure that 
the bounds of the classical region are correctly specified throughout the simulations 
- in which case the theory is referred to as the Projected Gross-Pitacvskii Equation 
(or PGPE). Starting from appropriately selected (but largely irrelevant) initial non- 
equilibrium conditions, one studies the system evolution to equilibrium; this approach 
is non-perturbative, and has even been used by some researchers at /near the critical 
region. 

(ii) The Truncated Wigner Approximation: This approach additionally includes 
quantum fluctuations, but these are only introduced in the initial conditions of the 
simulation. The subsequent evolution is again based on the (P)GPE, which here 
arises however only as an approximate equation for the full quantum evolution of the 
Wigner quasiprobability distribution function. As a result, such evolution actually 
leads to spurious damping, unless the study is limited to relatively short times or 
low temperatures. Contrary to the classical field methods, this approach is good for 
studying circumstances where quantum effects dominate the behaviour of the system, 
e.g. in optical lattices. 

(iii) The Stochastic Gross-Pitaevskii Equation (SGPE): This approach combines 
all the 'essential' elements of the successful previous theories (self-consistent Gross- 
Pitaevskii-Boltzmann or 'ZNG', classical field and Truncated Wigner) into a 'single- 
package', whose full potential has yet to be numerically unleashed. The SGPE 
resembles the ordinary GPE, but with the addition of a dissipativc term (as 
appropriate for perturbative mean field theories) and a dynamical noise term. In 
brief, the SGPE describes the dynamical coupling of the low-lying modes of a trapped 
gas to the 'reservoir' of higher-lying (thermal) modes, including both coherent and 
incoherent mechanisms. In principle, the higher-lying modes also evolve according to 
the quantum Boltzmann equation which further accounts for particle redistribution 
within the thermal cloud, as well as particle-exchanging collisions. 

8.4-1- Comparison between above Approaches: The stochastic formulation has the 

formal advantage over the various mean field theories that it is an explicitly number- 
conserving theory, which additionally includes a stochastic term. This theory 
essentially reduces to a description in terms of a dissipative Gross-Pitaevskii equation, 
in which the condensate energy is treated as a differential operator, c;oupled to a 
Quantum Boltzmann equation. Removing the operator nature of this energy, and 

tt Note that our discussion omits two other powerful theories (Positive-P Representation and 2- 
Path-Irreducible Effective Action) which arc currently gaining ground against some of the above 
approaches, and whose main added benefit is that they can handle strongly-correlated regimes that 
this Tutorial was not concerned with. 
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imposing a 'symmetry-breaking' approximation leads to a semi-classical theory which 
Zaremba, Nikuni and GrifRn have termed the 'ZNG' theory. Note that the latter 
theory introduces the additional simplification that interaction energies and chemical 
potentials are not treated fully self-consistently, but only to first order in the effective 
interaction strength. However, unhke that latter theory, instead of dealing with 
a coherent condensate, the SGPE provides a different physical interpretation, with 
$(r,t) referring to the full matter- wave field; this additionally includes both thermal 
and quantum fluctuations, and thus encompasses both the classical field theory and 
the Truncated Wigner method. In fact, the early formulation of the POPE was as the 
limiting case of an equation for the low-lying modes of the system when the coupling 
to a reservoir of higher- lying modes normally present is artificially removed. Inclusion 
of this coupling would essentially lead to the SGPE. 

The role of the coupling of the low-lying modes of the system described by the 
SGPE is to ensure that the system relaxes to the correct thermal equilibrium. Due to 
its coupling to the reservoir of higher-lying modes, the SGPE does not suffer from the 
low-temperature limitations of Truncated Wigner which is prone to spurious damping 
(heating) as a result of classical thermalization during the simulations. Although by 
its mathematical construction the SGPE formalism, coupled to a Quantum Boltzmann 
Equation, is in our opinion the most advanced formalism to date, its full potential has 
not yet been explored. In particular, in terms of its numerical implementation, most 
simulations to date have been performed (i) in the classical approximation, i.e. the 
thermal modes are not treated by the full Bose-Einstein distribution function, and 
quantum effects are discarded, and (ii) in coupling to a static heat bath, rather than 
to a dynamical thermal cloud. It is anticipated that both of these difficulties will be 
gradually overcome in the future. 

8.5. The Quest for an 'Optimal' Theory 

So, where does that leave us regarding a choice for an 'optimal' theory? We hope that 
the preceeding discussion has convinced the reader that the answer to this question 

actually depends on the details of the particular problem under consideration. While 
the most advanced theories can also be used to study much simpler problems, clearly 
one normally seeks the simplest (and computationally fastest) approach that will 
describe each particular case; we have thus attempted an approximate - but neither 
unique, nor free from controversy - classification along those lines in Fig. 16. In brief: 

• If dynamics are not an issue, then the self-consistent Hartree-Fock theory is a 
very good starting point, with Hartree-Fock-Bogoliubov-Popov or generalized 
Hartree-Fock-Bogoliubov perhaps constituting a 'safer bet'. However, if one is 
investigating low-dimensional systems at equilibrium, one should fully account 
for phase fluctuations by means of the appropriately modified mean field theory. 

• If one is interested in coupled condensate-thermal cloud dynamics (beyond the 

linear response limit), then one should in first instance resort to a 'two-component' 
formulation in terms of a dissipative Gross-Pitaevskii equation, coupled to a 
Quantum Boltzmann equation - sometimes referred to as the 'ZNG' theory. 

• If one is interested in studying regimes with large phase fluctuations, e.g. 
effectively low-dimensional geometries, spontaneous processes and related 
phenomena, or even the physics of the phase transition itself, then one should 
resort either to classical fleld or to stochastic methods. An added feature of these 
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approaches, which are computationally more demanding, is that they do not make 
an 'unphysical' separation of the atomic gas into condensate and thermal cloud 
(which is nonetheless desired by many researchers for ease of visualization), but 
rather they produce images reminiscent of experimental density profiles. The 
task of extracting information from such theories is thus somewhat more involved 
than for other types of theories, as one typically needs to perform a similar data 
analysis as would be done for analyzing an experiment (e.g. bimodal fits, averages 
over different realizations). 

In light of our preceeding analysis, our personal bias is to recommend a description 
based on the Stochastic Gross-Pitaevskii Equation for the predominantly coherent 
modes of the system, coupled to a Quantum Boltzmann Equation for the higher-lying 
modes of the system, a view generally shared by numerous - but certainly not all 
- researchers. Such a 'complete theory' has not yet been numerically implemented, 
although all essential features have been individually tested. Here we note that (i) the 
Stochastic Gross-Pitaevskii equation has been numerically simulated in coupling to a 
static heat bath by Stoof and collaborators (including one of us - NPP) and by Davis- 
Gardiner and collaborators, with all such treatments currently limited to one or two 
dimensions; (ii) moreover, a dissipative Gross-Pitaevskii equation, which is closely 
related to the stochastic Gross-Pitaevskii equation in the absence of the stochastic 
noise term, has been numerically solved (also by the present authors) in full coupling 
to a dynamical thermal cloud in three dimensions beyond the ergodic approximation. 
We believe that the implementation of such a theoretical framework constitutes the 
next significant development in the modelling of finite temperature trapped Bose gases. 

The issue of non-equilibrium dynamics of weakly-interacting Bose gases remains 
a fascinating topic of research, with a diverse range of theories covering different 
elements of relevance to current experiments. More challenging experiments are bound 
to be undertaken in the near future, particularly at/near the transition point, and the 
comparison of existing / improved theoretical models to such experiments over the 
course of time will determine which (if any) of these theories will live up to their 
expectation and whether one day one of these theories may be 'universally' recognized 
as the 'standard' theory for such systems. 
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Appendix A. The JILA Kinetic Theory 

The extended beyond-HFB second order perturbation theory discussed in Sec. 4.4.1 
was originally developed by Walser, Williams, Cooper and Holland in 1999 [110, 112]. 
Their treatment is based on the construction of a generalised kinetic theory for 
a 'coarse-grained' many-particle density operator which depends only on a few 
carefully selected 'master' variables; these were identified as the condensate mean 
field and gaussian fiuctuations around it, i.e. normal and anomalous averages, which 
are precisely the HFB variables. Like the treatment presented in Sec. 4.4.1, their 
derivation relies on a separation of timescale argument, under the assumption that 
the duration of a typical coUisional event is much smaller than the inverse collision 
rate; this enables the development of a systematic perturbative expansion. Below we 
briefiy review this formalism, highlighting the fact that all expressions are formulated 
in terms of a single-particle energy basis. The notation presented here explicitly 
maintains spatial dependence and preserves geometrical transformation properties. 

The Bose field operator ^{v,t) is expanded in a complete single-particle basis 
|1) via ^{T,t) = X^i(T'|l)ai- The binary interaction hamiltonian can be expressed in 
symbolic compact notation as 

H = {l\ho - m|2)*I*2 + yi234*I*i*3*4 , (A.l) 

where an implicit summation over repeated indices is assumed. Here (l|/io — 
/x|2) = /dr(l|r)[/io(r) - M](r|2), and ^1234 = (<?/2) / dr(l|r)(2|r)(r|3)(r|4). The 
latter expression has made use of the pseudopotential approximation, with resulting 
expressions being symmetrised, i.e. P1234 = V1243 = V2134 = ^2143- We introduce the 
slowly-varying 'master' variables describing the system as follows [110, 112]: 
The condensate is defined by 

|<A) = </.i|l) = (*) (^t) = (<^| = ^*(i| (A.2) 



CONTENTS 



105 



The normal average is defined by 

/ = (*t^,) = r + f= |<^)(</,| + /i2|l)(2| (A.3) 

where /° the coherent and / the incoherent contribution (the latter quantity is the 
analogue of p defined by Eq. (93) and used in Sec. 4.4.1). This notation can be 
understood by expressing the normal average in position representation via 

/(ri,r2) = (ri|/|r2) = {^{r^Mri)) = 0*(r2)</'(ri) + {SHr2)S{r,)) 

= (riirirs) + (ri|/>2) = ^(ri.rs) + /(n.rs) . (A.4) 

Similarly, the anomalous average corresponding to pair correlations is defined by 

m= (M) = 771^= + ^,= + mi2|l)|2) (A.5) 

where m'^ the coherent and m the incoherent contribution (analogue of k of Eq. (93)). 
In position representation, we find 

m(ri,r2) = (ri|(r2|m| = (^(r2)*(ri)) = 0(r2)0(ri) + (5(r2)<5(ri)) 

= (ri|(r2|m'' + (ri|(r2|?n = m''(ri,r2) +m(ri,r2) . (A.6) 

As in Sec. 4.3.2, we express the coherent and fiuctuation parts of the system (/, 
to) in generalised matrix notation via (c.f. expressions of Rc and Rnc in Eq. (110)) 



r + i 

The corresponding generalized hamiltonians for the condensate 11 (c.f. He) and the 
non-condensate S (c.f. Hj^c of Eq. (Ill)) are defined by 

where we have introduced 

Un = ho- n+Ufc +2Uf IIa = U^ 

T,N = ho-iJ. + 2C//C + 2Uf Y.A = Um (A.9) 

with N and A labelling normal and anomalous averages. The above equations make 
use of the shorthand notation (c.f. expressions for h and A of Eqs. (94) and (95)) 

Uf = 2Fi2'3'4'/3'2'|l)(4'| = 2yi2'3'4' (/3^.2' + h'2) |1)(4'| 

Um = 2Fi2'3'4'my4'|l)|2') = 2yi2'3'4' {ml,^, + TO3'4') |1)|2')| • (A.IO) 

The full second order dynamical theory which includes anomalous averages can 
thus be cast in the form (setting h= \) 

j^G< = -i^G< + {T<G> - V>G<) + h.c. . (A.ll) 

(Note the slight differences in notation between [112] and [113, 123].) In their simplest 

first orrfer limit, where only the first term (and its conjugate) is maintained, these give 
rise to the time-dependent HFB equations of Eq. (112). The remaining contributions 
appearing in Eq. (A.ll) include all second order coUisional integrals (expressed in 
terms of single-particle operators) . In particular, we have defined 



N J \ 



(A.12) 



the vectors, defined by cti 
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with the time-reversed matrices generated from these via = —a-iY^*ai and 
r> = — (Tir^*(Ti, where Gx the usual Pauli matrix exchanging the two components of 

1 

1 

Considering the condensate evohition first, the forward and backward transition 
rates {Y^, Y^, Y^ and Y^) describe the scattering of non-condensate particles into 
and out of the condensate; they are given by the following scattering processes F: 

Yn — ^ff{f+l) + '^^fmm* — ^fhmfh'- + fm{f+\) (A. 13) 

with the 'out' rates denoted by the superscript '>' defined by the corresponding 
expressions, obtained directly from the above by replacing (/+!) by / (and vice versa) 
in all expressions. This highlights the fact that only normal fluctuations / become 
bosonically enhanced, with the condensate and anomalous fields not exhibiting such 
an effect. 

The rates entering the corresponding expressions for the non-condensate axe 

= ^fhf+i) + ^/7'=(/+i) + ^fff 



+ 2 



2 



(A. 14) 



Note that the first contribution, r^^^j^j^-j, to explicitly contains the quantity 

f = f + f (unlike the remaining contributions which axe expressed in terms of /"^ or 
/). The corresponding '>' rates are obtained by replacing / by (/ + 1) and vice versa 
(and correspondingly /■*-»(/+ 1))- The collisional processes introduced above are 
defined by 

I"/// ~ 8Vi2'3'4'Vl"2"3"4"/3'l"/4'2"/4"2'|l)(3"| 
T/m/ = 8Vi2'3'4'Vi"2"3"4"/3'l"m4'3"/4"2'|l)|2") 
T/mm* = 8Fi2'3'4'Vl"2"3"4"/3'l"?Tl4'3"'Tl2"2'|l)(4"| 

rmmm* = 8Fi2'3'4' Vl"2"3"4"m3'4"m4'3"m2//2' |1) 1 1") , (A.15) 

where we have introduced the notation 

1 



14"2"3"4" = Vi"2"3"4" 



TrS{As") + iV 



Ae" 



where Ae" = e" + £3 ~ £3 ~ £4 denotes the difference in the single-particle energies 
between incoming and outgoing states. 

To gain some insight into the above lengthy expressions, we now highlight those 
contributions to the collisional integrals which involve only 'normal averages' / or /° 
(i.e. we ignore anomalous averages). We thus findff 

= ■ ■ ■ + (^/7(/+i) ~ (/+!)/) ^ {k.lQ) 

cl f 

= • • • + r//(/+i)(.f + 1) - F(^_^i)(j^i)j/ . (A.17) 

tfFor simplicity, we ignore here the second and third contributions to F^. 
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These equations contain the multi-mode generaUsations of the condensate and non- 
condcnsatc evolution given within our 'toy model' of Sec. 4.4.1. In particular: 

(i) the expression for the condensate mean field </> corresponds to Eq. (122). 

(ii) Upon noting that / = + Eq. (A. 17) can be re-expressed as 

+ r/e/(/+i)(/ + l)-r,.(;+i);/. (A.18) 

These two contributions correspond respectively to Eqs. (123) and (124), i.e. the 
collisional integrals of the Quantum Boltzmann Equation describing both population 

transfer collisions between condensate and thermal atoms, and collisions between non- 
condensate atoms which lead to thermal (incoherent) population redistribution. 

Appendix B. Method of Non-commutative Cumulcints 

The full system evolution can be derived from the fully non-local Hamiltonian of Eq. 
(15) via the study of suitably truncated equations of motion for averages of the field 
operator \l/(r,f), and its products: 

zfi|(*(r))(i) = ([*(r),if])(t) 



zfi^(#t (r) ■■■^{vm^i [¥{v) ■ ■ ■ #(r), ij] )(i) 

A possible consistent truncation scheme was discussed by one of us (NPP) in [97], 
whereby all averages of products of up to three fiuctuation operators were maintained, 
with averages of products of more than three operators decorrelated into products 
of averages of up to cubic terms. Such equations are to be solved self-consistently, 
with the requirement that the jkill interatomic potential is maintained. (It would 
be incorrect to use a pseudopotential here since it is precisely the role of some of 
these correlations to upgrade the exact interatomic potential to an effective one). 
This Tutorial has focused on the Markovian limit of such equations, where collisional 
processes are assumed to occur much faster than the typical evolution timescales 
of mean fields, thus rendering the state of the system independent of any initial 
correlations that may have arisen at some (distant) time in the past, i.e. there are 
no 'memory' effects. While this may be true for a large range of experiments with 
ultracold atoms, this condition is clearly violated in the experimentally-relevant case 
of attractive condensates [19, 20] or in the vicinity of Feshbach resonances [21], 
where a virtual bound state is supported close to the dissociation threshold of a 
two-body potential, both of which have been realized in recent experiments. (Note 
however that a discussion on non-Markovian effects was given in the context of the 
formalism of Sec. 4.4 in [115]). The method of non-commuting cumulants [297, 298] 
provides a mathematically well-defined procedure for incorporating such effects within 
a finite temperature scheme based on a suitable truncation of the hierarchy of coupled 
equations of motion. 

Mathematically, the cumulants can be defined as functional derivatives of a 
generating functional via [299] 

(*t(r„)...^-(ri))e 

S S 



(5J(r„) W*(ri) 



ln(e/'^'--^*W*('-'*)+^W*'(--'*))|j=j.=o , (B.l) 
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in the limit when the externally-added symmetry breaking field J(r) goes to zero. The 
cumulants may be obtained recursively (for general operators Oi, O2, O3, • • •) via 

(Ol) = (Ol)e 

(OlOa) =(Oi02)c+(Ol)c(02)e 
(616263) = (Oi0203)c + (Ol)c(02>c(03)c 

+ (Oi)e(0203)c + (02)c(Ol03)c + {O^i) c{Ol02) c (B.2) 

Cumulants of order higher than two provide a measure of how far the system is from 
the ideal Bose gas in thermal equilibrium (where such higher order correlations vanish). 
An approximate method to truncate the hierarchy of coupled equations of motion for 
the non-cummutative cumulants was proposed by Fricke [297]: Expanding the field 
operators in terms of single-mode operators one can decompose products of multiple 
operators into sums of pairwise contractions by means of Wick's theorem. Writing 
down the time-dependent equations for the correlation functions up to a given order of 
cumulants yields a closed system of coupled equations. The simplest possible scenario 
would be to include only first order cumulants (n = 1), thus obtaining the following 
evolution for the condensate mean field defined via ^'c(r,t) = (^(r))c(t): 



ho + J dr'F(r-r')|M/c(r',i)pM/c(r,t) 



(B.3) 



This equation should not be misinterpreted as the GPE of Eq. (29), as it is explicitly 
expressed here in terms of the actual interatomic potential V{r — r'). The required 
upgrade to an effective interaction (by the inclusion of multiple scattering) can be 
performed by the following self- consistent scheme: For a cumulant of order n, one 
must include the 'free dynamics' of normal-ordered cumulants of orders (n -|- 1) and 
(n + 2), i.e. their evolution for which cumulants of order (n + 3) and (n + 4) are 
ignored within their respective equations of motion. We follow Kohler and Burnett 
in explaining this in more detail [298]. Assuming we are mainly interested in the 
condensate dynamics (i.e. n = 1), consistent implementation of this scheme requires 
a systematic consideration of the cumulants of up to order (n + 2) = 3. The exact 
dynamics of the condensate becomes [298] 

d 

ih—^c{r,t) = /io(r)*c(r,t) 

+ J dr'V{r-r')^l{r',t) [*e(r, t)*c(r', *) + $c(r, r', t)] 
+ j drV(r-r') [^,{r' ,t)r,{r,r' ,t) + ^,{r,t)r,{r' y ,t)] 

+ J rfr'l/(r-r')(*^(r')*(r')*(r))c(t) , (B.4) 

where we have defined the second order cumulants $c(ri,r2,t) = (^(r2)^'(ri))c(t) 
(pair function), and rc(ri,r2,i) = (*I'^(r2)^'(ri))c(t) (non-condensate one-body 
density matrix) (at the level n + 1 = 2). Note that the third order cumulant 
{^^r')^ {r')^ (r)) c{t) also appears in this expression, consistent with Eq. (117). 

To proceed further, one should formulate equations of motion for these cumulants 
to the desired level of truncation. We find [298] 

i/i^$e(r, r', t) = [ho{r) + ho{r') + V{r - r')] $c(r, r', t) 
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+ y(r-r')*c(r,t)*c(r',i) , 



zfi-r,(r,r',t) = \ ha{v)-ho{v')\ T^{v,v',t) , 

= 'ho{r) + Ao(r') + V{r - r') - ho{r")] (*t(r")*(r')*(r))e(i) 

+ V{r - r') [*e(r', i)r(r, r", t) + *e(r, i)r(r', r", t)] . (B.5) 

(In obtaining the corresponding equation of motion for {r")'^{r')'^{r))c.{t) we 
have ignored products of normal-ordered cumulants containing n + 4 = 5 Bose field 
operators.) Solving these equations formally by means of two-body Green's functions 
(in the basis of trap eigcnstatcs) generates an cxac;t crloscd non-Markovian nonlinear 
Schrodinger equation for the condensate wavefunction, in which the interaction term 
has been formally renormalised to the two-body T-matrix - see Ref. [298] for detailed 
expressions and further details. Note that the usual GPE of Eq. (29) arises as a 
limiting case of this equation upon imposing the Markov approximation and using 
the contact potential approximation for the two-body T-matrix. The above procedure 
generalises straightforwardly for larger n: for example, for n = 2 one recovers a 
generalised non-local non-Markovian form of the HFB equations. 

The method of non-commutative cumulants has been successfully applied to 
diverse problems including the scattering of two condensates [298] and the atomic- 
molecular BEG problem in the vicinity of Feshbach resonances [21] - for a recent review 
see [25]. An important advantage of this approach is that it can incorporate three- 
body interactions by means of non-local pairwise potentials, thus yielding generalised 
three-body T-matrices satisfying Fadeev equations [300]. 



Appendix C. The Gctrdiner-Zoller Kinetic Theory 

Appendix C.l. The Condensate Band Master Equation 

The full master equation for the evolution of the condensate hand density operator 
developed by Gardiner and ZoUer is given by [214, 130, 215, 216, 217] 

= pH^'" + pGrowth ^ ^Scatt _ 

The first contribution arises entirely from internal condensate band dynamics, via 
^Ham _ _(j/;j)[if^^ p^]^ where is the condensate band hamiltonian which also 
includes the effect of the average non-condensate density on the condensate. The 
remaining two terms arise from the interaction between the condensate and the non- 
condensate band, and are obtained by consideration of the interaction hamiltonian 
i^int. In presenting the explicit form of these contributions, we use for simplicity 
the notation — (j) + 5, although it should be understood that these operators are 
actually the corresponding number-conserving operators discussed in Sec. 5.2.2, which 
have been already implicitly projected onto orthogonal subspaces. To write down the 
full master equation, we additionally define the centre of mass and relative coordinates 
ro = {y + v')/2 and f = (r — r'), where r and r' the coordinates of the colliding atoms. 

Like in our beyond-HFB mean field discussion (Sec. 4.4), the growth terms 
corresponding to particle transfer between the two bands arise from the contributions 
to the system hamiltonian which involve three non-condensate operators 5"^ 56 (i.e. from 



CONTENTS 



110 



the corresponding number-conserving generalisation oiH^), which yield a contribution 
of the form 



•Growth 



= j dvoj dv {(^1 - A2) + (Bi - B2)} 



(C.2) 



Each of the four contributions appearing in the integrand are given in terms of a 
suitable commutator of a product involving the condensate band density matrix pc, 
a creation (annihilation) operator 5^^\ and a matrix or G^") corresponding to 
scattering rates into and out of the condensate band. In particular 



^1 = 



G(-)(Le)0|ro 

PeG(-)(-Lc) 

GW(-ic)^^ fro- 

PcG(+)(Xe) 



ro - 



r 

ro + 
ro + 



'■°-2 



r 

^°+2 



(C.3) 



Here G(+) and G^") depend both on ro and r (suppressed in above expressions) and 
on energy, with defined by Lc(i){v) = [0(r) , H^. They are given by integrals of the 
form 



G(±)(ro,f,a;) 



X e 



(27r)8fi2 

»(k2+k; 



(ik2 / (iks / d\s.4^5{uj2 + ui^ — uo^ — lo) 



'".^^(ro,ki) 



Here .F+(ro,ki) = FiF2(i^3 + 1) and J"-(ro,ki) 
the one-particle Wigner function introduced by 



(fl + l)(F2 



(C.4) 

1)-F3, and Fj denotes 



1 



(27r) 



dkF(ro,k)e-'(''"'°+'^*\ (C.5) 



where to — hk'^/2m + K!xt(ro) and the range of integration has been restricted to 
the non-condensate band. If the non-condensate band is in thermal equilibrium, 
the 'in' and 'out' factors are related (as discussed in Sec. 6.2.4) by G'-+-'(w) = 
g-/3(Ba;-M)(^(-)(-^,^ Note the strong similarity of the two terms of Eq. (C.4) to the 
contributions (corresponding to 'in' and 'out' rates) contained in iR{r, t) of Eq. (153) 
and 5]^(r) of Eq. (199). 

The scattering terms describe collisions of atoms in different bands, which lead 
to population redistribution within each hand, but cause no net transfer between 
bands. They arise from contributions to the system hamiltonian which contain one 
non-condensate creation (5^) and one non-condensate annihilation [S) operator (and 
also a product of condensate creation and annihilation operators) - i.e. effectively 
from the 'resonant' part of the non-condensate form of the number-conserving H2 
contribution. Their evolution takes the form 



• Scatt 



j dvo j df {Ci 



-C2} , 



(C.6) 



where 



ro 



ro 



ro 
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X 5{u}2 - T w)Fi (F2 + 1) . (C.7) 
Again, at thermal equilibrium M(ro,f,w) = e~^'^ M{vo,v, —lo). 

Appendix C.2. The Gardiner- Davis Stochastic Gross-Pitaevskii Equation 

Here we give explicit expression for the terms appearing in the local form of the 
Gardiner-Davis stochastic GPE of Eq. (214) [216]. The growth term is given by 
[215, 216, 217, 241] 

Kcir, t) = pG{r) (/x - HL,) a(r, t) (C.8) 

where G(r) = / drG'^+\r,r,uj = 0) with G'(+) defined in Eq. (C.4) , and is 
the condensate eigenvalue obtained from the projected finite temperature GPE via 
Lca(r) = T'c + 2,gn7vc(r) + .9|a(r)p| Q:(r). Here Vc is a projector onto the 

subspace orthogonal to the trap basis, defined via Pc{^,r') = 1 — X)„ (r)i^(r')' 
where Yn{r) are a complete set of trap eigenfunctions and the summation is restricted 
to modes of the non-condensate band. In the above expression, we have also defined 
the average non-condensate band density via njvc(i") = T^^nc 
corresponding scattering term, KM{^,t) is defined by [216] 

KM{r,t) = -^M(r)a(r) {a* {r,t)LMr,t) - a{r,t) [L,a{r,t)]*} ,(C.9) 

where the parameter M is an approximate local form of the quantity M(ro,f, ±0;) of 
Eq. (C.7). The correlations of the noise appearing in Eq. (214) are 

dWair, t)dWG{r', t) = 2G(r),5(r - v')dt 

dWG{v,t)dWG{v',t) = dW^{v,t)dW^{v',t) = 

rfWM(r, t)dWM{r', t) = 2M(r)(5(r - r')dt . (C.IO) 

We now demonstrate the equivalence between the stochastic equations of 
Gardiner, Anglin, Fudge and Davis (Eq. (214)) and that of Stoof (Eq. (200)). For 
this purpose, we firstly ignore the scattering contributions (which are cumbersome to 
deal with) and recast Eq. (214) in the simpler form 

1 

da{r,t) = — -VcLGpa{r,t)dt 

- [-if3hG{v)] Ve {Lgp - m) a(r, t)dt + dWG{v,t). (C.ll) 

Firstly, we eliminate the 'trivial' time-dependence by re-expressing the equation in 
terms of a{v,t) = e'^*/''a(r, f) to obtain 

da{r, t)= - ;^ {[1 - i/?fiG(r)] Vc {Lgp - m)} a{v, t)dt 

+ dWc{v,t), (C.12) 

where Lgp is the eigenvalue of the ordinary Gross-Pitaevskii equation, defined by 
Lgpa(r) = {ho + g\a{Y)\^)a{Y) and dWG{Y,t) = e^^'*/^dWG{Y,t). Upon noting that 
5^if(r) = — 4iG(r) we see that this equation is identical to Eq. (200). It only remains 
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to prove that the noise sources are the same; the noise contribution corresponding to 

the differential equation ihda{r,t)/dt becomes ihdWG{'r,t)/dt; using the expressions 
of Eq. (C.IO), its correlations thus become 

{(^-iMWG*{r,t)) (iMWG{r',t)))=h^dW^{r,t)dWG{r',t)) 

= 2!fG{r)Sc{r,r') ^ '-tfi:'' {r)5{r - r') , (C.13) 

demonstrating the direct equivalence between the stochastic equations of Stoof and 
Gardiner and co-workers in the Umit where the scattering processes can be ignored. 
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